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FOREWORD 


This  report  is  a  slightly  revised  version  of  the  dissertation  submitted  in  partial 
fulfillment  of  the  requirement  for  the  Doctor  of  Philosophy  degree  at  the  University  of 
Washington  in  December  1994.  Dr.  James  C.  Luby,  Principal  Engineer  at  the  Applied 
Physics  Laboratory  and  Research  Assistant  Professor  of  Electrical  Engineering,  was  the 
chairman  of  the  supervisory  committee.  Financial  support  for  this  research  was  provided 
largely  by  the  APL  Fellowship  Program,  with  supplemental  funding  by  NRL  Code  7105 
(Dr.  E.  R.  Franchi)  and  by  the  Environmental  Calibration  by  Acoustic  Holography 
Program  sponsored  by  the  Office  of  Naval  Research,  Code  3210A. 
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ABSTRACT 


In  some  underwater  acoustic  applications,  it  is  desirable  to  predict  what  an 
acoustic  field  will  be  at  some  distance  away  from  a  source.  Numerical  modeling  based  on 
environmental  parameters  (e.g.,  sound  speed,  bathymetry,  sediment  properties)  is  one 
possibility,  but  the  amount  of  environmental  information  necessary  for  accurate  long- 
range  modeling  is  often  beyond  what  is  realistically  measurable.  Actual  acoustic 
measurements  are  another  possibility,  but  our  knowledge  of  the  fields  is  then  limited  to  the 
points  in  the  ocean  where  the  sources  and  receivers  are  located.  We  present  a  hybrid 
method  for  predicting  acoustic  fields  that  takes  limited  acoustic  measurements  and 
extrapolates  them  over  short  ranges  using  modeled  fields.  The  measured  and  modeled 
fields  are  combined  in  an  integral  formulation  via  a  specialization  of  Huygens’  principle. 
The  basic  theoretical  formulation  is  stated  and  analyzed  using  the  theory  of  normal  modes. 
The  formulation  leads  directly  to  modal-dependent  (grazing  angle  dependent)  “obliquity” 
factors,  which  are  canceled  to  first  order  by  averaging  with  a  range  derivative-based 
extrapolation.  Sensitivity  of  the  algorithm  to  the  higher  order  components  of  the  obliquity 
factors  is  studied.  As  with  other  model-based  algorithms,  environmental  mismatch  can 
degrade  the  results.  The  algorithm’s  sensitivity  to  mismatch  in  water  column  sound  speed 
and  sediment  sound  speed,  attenuation,  and  density  are  studied  via  a  perturbational 
approach.  Since  the  algorithm  explicitly  contains  an  integral  and  the  measured  field  is 
collected  at  discrete  locations,  numerical  quadrature  techniques  are  studied.  It  is  shown 
how  undersampling  can  degrade  the  quality  of  the  extrapolated  fields,  and  how  the  spatial 
sampling  required  for  reliable  extrapolation  can  be  obtained.  A  modal  decomposition- 
based  integration  method  is  then  shown,  where  the  number  of  reference  elements  required 
is  equal  to  the  number  of  modes  that  contribute  significantly  to  the  fields.  A  modified 
version  of  this  formulation  can  exactly  cancel  the  obliquity  factors. 
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Chapter  1:  Introduction 


1.1.  Introduction  to  the  Dissertation  and  to  Chapter  1 

It  is  well  known  that  electromagnetic  waves  do  not  propagate  well  through  the 
ocean  due  to  the  strong  attenuating  effects  of  seawater  on  this  kind  of  energy.  On  the  other 
hand,  under  the  proper  conditions,  acoustic  waves  can  effectively  propagate  over 
extremely  long  distances  (on  the  order  of  thousands  of  kilometers).  Even  over  relatively 
short  ranges,  the  use  of  acoustic  energy  is  preferred  for  a  variety  of  underwater 
applications  such  as  communications,  active/passive  target  detection  and  localization, 
imaging  of  submerged  objects,  etc. 

In  some  underwater  acoustic  applications,  it  is  necessary  to  predict  what  an 
acoustic  field  will  be  at  some  point  distant  from  an  acoustic  source.  This  is  usually 
accomplished  by  first  sampling  the  environmental  properties  of  the  ocean  that  affect  the 
propagation  between  the  source  and  the  distant  point.  Examples  of  these  properties 
include  (but  are  not  limited  to)  the  sound  speed  in  the  water  column,  the  bathymetry  in  the 
region  of  interest,  and  the  sound  speed,  density,  and  attenuation  in  the  ocean  bottom.  For 
very  simple  situations,  analytic  solutions  of  the  acoustic  field  are  possible.  For  most 
scenarios  the  sampled  environmental  properties  are  used  as  input  to  one  of  several  types  of 
numerically-based  methods  for  modeling  underwater  acoustic  fields.  Two  such  methods 
currently  implemented  in  computer  software  are  parabolic  equation  (PE)  approximations 
[Collins  and  Westwood,  Collins]  and  normal  mode-based  solutions  [Porter  and  Reiss, 
Porter  (1992)]. 

In  practice,  however,  the  amount  of  environmental  information  required  for 
accurate  predictions  of  the  acoustic  field  over  any  significant  range  may  be  beyond  what  is 
realistically  measurable.  For  example,  the  bathymetry  and  the  sound  speed  profile  in  the 
water  column  may  vary  too  rapidly  with  range  and/or  time  to  allow  adequate  sampling. 
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Ocean  bottom  parameters  (which  are  typically  difficult  to  measure  in  any  case)  may  also 
vary  unpredictably  with  range.  Surface  and  bottom  roughness  can  also  contribute  to 
modeling  difficulties. 

Another  possible  method  for  determining  acoustic  fields  distant  from  a  source  is  to 
actually  measure  them  using  hydrophone  arrays.  Although  this  method  would  provide 
exact  results  for  any  kind  of  environmental  situation,  its  use  is  limited  by  the  fact  that  the 
fields  would  only  be  known  for  the  given  source  and  receiver  locations. 

This  dissertation  will  describe  and  analyze  a  hybrid  method  for  predicting 
underwater  acoustic  fields  that  utilizes  both  limited  acoustic  measurements  and  modeled 
fields  (over  relatively  short  ranges)  that  are  based  on  local  environmental  information.  The 
two  types  of  information  are  gathered  along  arrays  of  transmitters  and/or  receivers 
deployed  vertically  in  the  ocean,  and  are  “combined”  using  a  specialization  of  Huygens’ 
principle.  The  measured  acoustic  data  can  be  thought  of  as  “calibrating”  a  region  of  ocean. 
The  modeled  fields  are  then  used  to  “extrapolate”  the  measurements. 

The  remainder  of  Chapter  1  will  contain  an  introduction  to  normal  mode  theory, 
which  is  used  extensively  throughout  the  remainder  of  this  work.  It  will  also  review  some 
previous  work  on  a  holographic  array  processing  (HAP)  algorithm  for  source  localization 
whose  mathematical  underpinnings  are  related  to  those  of  the  acoustic  field  extrapolation 
algorithm.  A  preview  to  the  remainder  of  the  dissertation  will  follow. 

1.2.  Normal  Mode  Theory 

This  section  is  intended  as  a  brief  introduction  to  normal  mode  theory  for 
underwater  acoustic  propagation.  While  some  of  the  “highlights”  of  the  mathematical 
derivations  will  be  presented,  mathematical  rigor  will  be  eschewed  in  favor  of  physical 
interpretations  of  the  resulting  equations.  This  material  will  follow  the  treatment  found  in 
[Jensen  et  a/.],  although  other  excellent  rigorous  mathematical  derivations  can  be  found  in 
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[Boyles]  and  [Brekhovskikh  and  Lysanov]. 

1.2.1.  Range-independent  Normal  Modes 

The  inhomogeneous  acoustic  wave  equation  for  pressure  can  be  derived  from 
basic  physical  laws  [Jensen  et  al]  and  stated  as 

2 

=  /(r,0,  (1-1) 

^  cW 


where  p  is  density,  c  is  the  speed  of  sound,  p  is  pressure,  and  /(r,  t)  is  the  source  term  as 
a  function  of  space  and  time.  Several  assumptions  are  now  made.  First  we  assume  that  we 
will  be  dealing  with  narrow-band  sources  of  the  type  e  This  justifies  Fourier 
transforming  the  wave  equation  into  the  frequency  domain  and,  therefore,  obtaining  a 
time-independent  solution.  Secondly,  we  assume  a  cylindrically  symmetric  medium  about 
the  source,  and,  therefore,  will  use  a  cylindrical  coordinate  system.  Lastly,  we  assume  that 
p  and  c  depend  only  on  depth,  i.e.,  the  medium  is  range-independent.  Taking  these 
assumptions  into  account,  we  can  write  the  so-called  Helmholtz  equation  as 


rdrL 


dp{r,  z) 
dr 


]  *  ^<4 


1  dp(r,  z)1  CO 


Lp(z)  dz 


+  ^p(r,z)  =  - 


c\z) 


8(r)5(z-z^) 
2nr 


(1.2) 


where  r  is  the  range  coordinate,  z  is  the  depth  coordinate,  is  the  source  depth,  and  0  is 
the  radian  frequency  of  the  source. 

In  order  to  solve  Eq.  1.2,  we  use  the  technique  of  separation  of  variables  [Boyles]. 
This  technique  assumes  a  solution  that  is  a  product  of  range-only-dependent  and  depth- 
only-dependent  factors,  i.e.. 


p(r,z)  =  (|)(r)\|/(z). 


(1.3) 
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By  substituting  this  into  the  homogeneous  Helmholtz  equation,  the  so-called  modal 
equation  can  be  written  as 


P(z)- 


1 


dzipiz)  dz  J 


=  0 


(1.4) 


with  boundary  conditions 


¥  JO)  =  0 


and 


dz 


=  0; 


z  =  D 


(1.5) 


2 

is  called  the  separation  constant. 

The  boundary  condition  at  the  surface  (z  =  0 )  implies  what  is  known  as  a  free  or 
pressure  release  surface.  It  is  also  known  as  the  Dirichlet  boundary  condition.  It  is 
characterized  by  the  fact  that  acoustic  pressure  vanishes  along  such  a  surface.  The  lower 
boundary  condition  (at  depth  z  =  D)  implies  what  is  known  as  a  rigid  surface.  It  is  also 
known  as  the  Neumann  boundary  condition.  It  is  characterized  by  the  fact  that  the 
derivative  of  acoustic  pressure  (or,  equivalently,  particle  velocity)  vanishes  along  this  type 
of  surface.  It  should  be  noted  that  while  the  ocean  surface  can  behave  much  like  a 
Dirichlet  boundary  under  calm  conditions,  the  Neumann  boundary  condition  is  often  a 
poor  approximation  of  the  ocean  bottom.  The  Neumann  condition  is  still  used  as  a 
mathematical  nicety  for  analytical  work,  however,  and  more  realistic  bottom  models  can 
easily  be  incorporated  into  numerical  modeling  algorithms. 

The  modal  equation  is  now  in  the  form  of  a  proper  Sturm-Liouville  eigenvalue 
problem  [Boyles].  Several  properties  [Jensen  el  al]  of  this  type  of  problem  will  aid  in  our 
understanding  of  material  to  come: 

( 1 )  the  modal  equation  has  an  infinite  number  of  solutions,  similar  in  nature  to  the 

modes  of  a  vibrating  string. 
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(2)  each  mode  is  characterized  by  an  eigenfunction,  or  mode  shape  function, 

(where  the  mth  mode  function  has  m  zeroes),  and  an  eigenvalue,  or  what  can 
be  interpreted  as  a  horizontal  propagation  constant  or  horizontal  wavenumber,  , 

(3)  the  eigenvalues  are  all  distinct,  are  ordered  such  that  ^  j  >  ^2  >  •  •  •  > 

shown  that  all  the  eigenvalues  are  less  than  where  is  the  lowest  sound 

speed  of  the  problem  [Jensen  et  al], 

(4)  the  eigenfunctions  form  an  orthonormal  set  relative  to  p(z) ,  i.e., 

- — - dz=  {  ^  ,  (1-6) 

J  p(z)  0,m^n 

0 

and  hence  the  terminology  “normal  modes,” 

(5)  the  eigenfunctions  form  a  complete  set,  so  that  the  pressure  can  be  expressed  as 
a  sum  of  normal  modes: 


P(r,z)  =  X  (1-7) 

m  =  1 

After  some  mathematical  manipulation,  and  using  the  above  properties,  it  can  be 
shown  that  [Jensen  et  al] 


«„'■).  (‘-8) 

where  Hq\x)  is  the  zero-order  Hankel  function  of  the  first  kind.  Using  the  asymptotic 
(large  argument,  x  »  1 )  approximation  of  the  Hankel  function,  i.e.. 


(1.9) 
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we  can  finally  express  the  pressure  as 


P 


(1.10) 


Note  that  this  equation  can  be  simply  thought  of  as  a  linear  combination  of  the  mode 
functions,  i.e., 


p{r,z)=  2  (1-11) 

m  -  1 

where  the  W^^’s  are  complex-valued  weights  that  are  functions  of  the  horizontal 
wavenumbers  and  of  the  mode  functions  at  the  transmit  depth,  density,  and  range. 

In  accordance  with  property  (2)  above,  each  mode  can  be  thought  of  as  having  a 
wavenumber 


where  |A:^|  =  for  all  m,  and  Ti^  can  be  considered  horizontal  and  vertical 
wavenumbers  respectively,  and  r  and  z  are  unit  vectors  in  range  and  depth  respectively. 
Each  mode  can  be  thought  of  as  describing  energy  that  propagates  at  an  approximate 
grazing  angle  (i.e.,  with  respect  to  the  horizontal) 


y  ~  acos 


(1.13) 


where  is  a  reference  wavenumber.  For  this  interpretation  to  make  physical 

sense,  should  be  equal  to  in  accordance  with  property  (3)  given  above.  Also  in 
accordance  with  property  (3)  is  the  fact  that  higher  order  modes  propagate  at  increasingly 
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Steeper  grazing  angles. 

1.2.2.  Weakly  Range-Dependent  (Adiabatic)  Normal  Modes 

The  previous  section  dealt  with  a  range-independent  medium.  When  the  medium 
varies  with  range,  a  phenomenon  called  mode  coupling  occurs.  This  means  that  as  the 
acoustic  energy  propagates  in  range,  energy  from  one  mode  can  be  transferred  into 
another  mode.  Variations  in  sound  speed  profiles  with  range,  sloping  bottoms,  internal 
waves,  surface  and  bottom  roughness,  and  other  factors  can  all  contribute  to  mode 
coupling. 

The  degree  of  mode  coupling  between  modes  m  and  n  at  a  given  range  r  is 
governed  by  the  relations  (Jensen  et  al.) 


mn 


fii- 

r 

J  rdr 

L  dr  J 

dz 


=  \-K - dz. 

J  or  p 


(1.14) 


The  so-called  adiabatic  approximation  assumes  that  these  quantities  are  negligible  and  the 
pressure  is  written 


p{r,z)  s 


1/2 


i 

SrcrJ  p(z ) 


(1.15) 


m  =  1 


The  mode  functions  and  horizontal  wavenumbers  at  any  given  range  are  determined  by 
using  the  reference  waveguide  method  [Brekhovskikh  and  Lysanov].  This  method  consists 
of  taking  the  sound  speed  and  density  profiles  and  the  depth  at  the  given  range  and 
constructing  a  hypothetical  range-independent  waveguide.  The  method  of  the  previous 
section  is  then  used  to  solve  for  the  modes  and  horizontal  wavenumbers  for  that  range. 
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Note  again  that  the  pressure  can  be  expressed  as  a  linear  sum  of  complex-weighted  mode 
functions  for  any  given  receive  range. 

It  can  be  seen  that  there  is  not  a  great  deal  of  difference  between  Eq.  1.15  and  its 
range-independent  counterpart  Eq.  1.10.  One  difference  is  the  accumulated  phase  in  the 
exponential.  The  integrated  phase  in  Eq.  1.15  averages  over  range  the  contribution  of  the 
range-dependent  horizontal  wavenumbers  to  the  total  phase.  Another  difference  is  the 
range  dependence  of  the  other  arguments  in  the  summand.  Note  that  Eq.  1.15  reverts  to 
Eq.  1.10  for  range-independent  mode  functions  and  horizontal  wavenumbers. 

1.2.3.  Stepwise  Coupled  Modes 

A  coupled  mode  solution  for  a  stepwise  range-dependent  environment  as  pictured 
in  Fig.  1.1  can  be  given  as  [Evans] 


Pfr,z)  =  +  (1.16) 

m 

where  j  denotes  the  range  segment  (;  =  1  is  to  the  left  of  r  j ),  m  is  the  mode  index,  the  A’s 
and  5’s  are  complicated  expressions  dependent  on  the  source  condition  and  mode 


Figure  1.1.  Schematic  for  stepwise  coupled  modes. 
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coupling  matrices,  and  the  HVs  and  H2's  are  ratios  of  Hankel  functions  given  by 


H, 


(1) 


= 


(^-  r) 


H, 


(2) 


ii-  r) 


(1.17) 


r^._  j  =  Tj  for  ;■  =  1  is  a  special  case.  Refer  to  [Evans]  for  a  more  detailed  explanation. 
Again,  note  that  for  any  given  range  segment  the  pressure  is  a  linear  combination  of  the 
reference  waveguide  mode  functions  for  that  range. 


1.3.  Holographic  Array  Processing  for  Source  Localization 

An  important  research  topic  in  recent  years  has  been  the  underwater  source 
localization  problem,  i.e.,  given  a  remote  measurement  of  an  acoustic  signal  generated  by 
an  unknown  source,  how  can  one  determine  the  location  of  the  source?  A  set  of  algorithms 
falling  under  the  general  heading  of  “matched-field  processing”  has  arisen  in  an  effort  to 
solve  this  problem.  The  algorithms  essentially  consist  of  comparing  a  measured  acoustic 
field  P  due  to  the  source  with  simulated  replica  fields  P^  that  are  numerically  calculated 
based  on  available  environmental  information.  The  replica  fields  are  calculated  for  a  grid 
of  hypothetical  source  locations,  and  the  P^  that  most  closely  “matches”  the  measured 
field  P  indicates  the  location  of  the  source. 

Consider  a  vertical  receiving  array  at  range  r  =  0  with  receivers  at  a  set  of  depths 

-iidt  .  , 

z  .  For  a  source  at  range  r,  and  depth  z,  emitting  a  time-harmomc  e  signal,  the 
measured  complex  pressure  at  receiver  n  is  denoted  P(r^,  z^,  z„)  •  Likewise,  the  replica 
field  for  a  hypothetical  source  position  (r^',  z/)  is  denoted  P^{rJ,  z^',  z„).  For  a  grid  of 
hypothetical  source  locations,  an  “ambiguity  function”  can  be  computed  as  [Hamson  and 
Heitmeyer] 
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A(r^,z^,r/,z/)  = 


SI/’I  X!'’J 


(1.18) 


where  the  *  represents  complex  conjugate.  The  ambiguity  function  will  be  one  for  the 
correctly  located  source  position,  and  less  than  one  for  incorrect  source  locations.  This 
assumes,  of  course,  that  the  modeled  field  corresponding  to  the  correct  source  location  is  a 
perfect  representation  of  the  true  field.  Other  algorithms  [Baggeroer  etai]  incorporate  the 
cross-spectral  matrix  for  the  array  into  the  algorithm,  but  Eq.  1.18  will  suffice  for  material 
presented  here. 

Calculation  of  the  replica  fields  z/,  z^)  is  often  the  most  troublesome  aspect 

of  matched-field  source  localization  algorithms  for  reasons  stated  in  Sec.  1.1.  Several 
studies  have  shown  how  insufficient  and/or  incorrect  knowledge  of  the  parameters  that 
affect  the  propagation  and  that  are  used  as  input  to  the  numerical  codes  can  degrade  the 
performance  of  matched-field  processors  (see  references  given  in  [Baggeroer  et  al.]). 
Typically,  as  the  range  over  which  one  needs  to  model  the  propagation  grows,  so  does  the 
error  that  is  introduced  into  the  calculated  replica  fields.  This  is  because  more 
environmental  variability  is  introduced  into  the  problem  as  range  increases,  meaning  that 
more  environmental  variability  must  be  taken  into  account  by  the  modeling  algorithm. 

It  was  observed  [Rouseff  (1989)]  that  one  could  use  variable  depth  reference 
sources  in  the  vicinity  of  the  unknown  source  to  partially  compensate  for  the  propagation 
path  between  the  unknown  source  and  a  receiving  array.  Figure  1.2  shows  the  general 
geometry  of  the  problem.  This  idea  was  then  formalized  [Mourad  et  al.  Porter  et  al 
(1992)]  into  a  source  localization  algorithm.  Using  adiabatic  normal  mode  theory  and  the 
notation  of  [Mourad  et  al],  we  can  write  the  pressure  field  incident  on  the  receiving  array 
due  to  the  unknown  source  as 
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fR 


Figure  1.2.  Geometry  for  source  localization. 


'’,(^..0)=  ^  Sv,(z,,0)v/Z,,r,)[4,.(/-,)J  '"'e  J"  (1.19) 

L  J  y 


and  that  due  to  the  reference  source  as 


(1.20) 


Here  a  constant  density  of  p  =  1  is  assumed,  and  a  slightly  different  scaling  factor 
(consistent  with  that  found  in  [Brekhovskikh  and  Lysanov])  is  used  in  place  of  the  one 
found  in  Section  1 .2.2. 

The  key  to  the  source  localization  algorithm  [Rouseff  (1989)]  is  that  we  can  form  a 
“filtered”  version  of  P,  by  calculating 


h{0) 

S  {z^,  r^;Zj^,  r^)  =  J  Pj^iz,  0)P,(z,  0)dz, 
0 


(1.21) 


where  h  (0)  is  the  depth  at  the  receiving  array.  Using  the  orthonormality  of  the  modes 
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allows  us  to  write 


5(z,,  r^)  = 


X 


(1.22) 


Due  to  the  form  of  the  eigenfunctions  and  the  accumulated  phase  (in  the  exponential)  in 
Eq.  1.22,  the  function  S{Zj,  r^)  can  be  interpreted  as  an  approximation  of  the  field 
due  to  the  source  at  (z^,  r^)  that  would  be  incident  on  an  array  located  at  the  range  of  the 
reference  sources.  We  can  say  that  by  using  the  depth-variable  reference  source  we  have 
created  a  “virtual”  receive  array  at  range  . 

The  holographic  field  S  now  allows  us  to  localize  the  source  using  a  matched-field 
approach  such  that  the  range  over  which  the  propagation  must  be  modeled  is  greatly 
reduced.  Specifically,  we  can  calculate  replica  fields  K{zj^,  rf^',z,  r')  at  the  location  in 
range  of  the  virtual  array  using  local  estimates  of  the  propagation  parameters  for 
hypothetical  source  locations  (z',  r')  .  Then,  as  an  analogue  to  Eq.  1.18,  we  calculate 


h{r^) 

Q(z',r')=  J  K{zj^,r^;z',r')S  iZj,r^;Z]^,rj^)dZ]^.  (1.23) 

0 


It  can  be  shown  [Mourad  et  al.]  that  the  Q  function  will  attain  a  maximum  value  as 
(z,  r')  -4  (z^,  r,)  . 

A  treatise  on  this  “holographic  array  processing”  (HAP)  algorithm  for  source 
localization  can  be  found  in  [Al-Kurd].  The  contribution  of  this  work  is  on  several  fronts. 
First,  it  contains  a  detailed  examination  of  the  focusing  condition.  Next  it  treats  the  case 
where  one  or  both  of  the  arrays  (receive  and/or  virtual)  are  truncated,  i.e.,  do  not  span  the 
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entire  water  column.  It  can  be  shown  that  truncation  of  the  virtual  array  has  much  more  of 
a  de-focusing  effect  on  the  algorithm  than  truncation  of  the  receive  array.  Finally,  the 
performance  of  the  algorithm  in  ambient  noise  is  examined,  and  specifically  how  it  is 
affected  by  the  SNR  of  the  reference  sources. 

1.4.  Preview  to  the  Dissertation 

Chapter  2  of  this  dissertation  will  outline  the  basic  theory  behind  the  acoustic  field 
extrapolation  algorithm.  It  will  start  by  discussing  a  possible  situation  in  which 
extrapolated  acoustic  fields  may  be  desired.  It  will  then  examine  some  theoretical  material 
from  the  field  of  optics,  namely  Huygens’  principle  and  the  Fresnel-Kirchhoff  diffraction 
formula,  which  will  serve  as  motivation  for  the  field  extrapolation  algorithm.  The  idea  of 
acoustically  “calibrating”  an  underwater  environment  will  be  briefly  discussed.  Assuming 
a  range-independent  medium,  the  extrapolation  algorithm  will  then  be  stated.  It  will  be 
shown  how  combining  fields  that  are  modeled  over  the  extrapolated  distance  with  the 
calibration  data  via  a  specialization  of  Huygens’  principle  can  produce  an  approximation 
to  the  actual  field  that  would  be  incident  on  a  receiving  array.  It  will  be  seen  how  modal- 
dependent  “obliquity  factors”  distort  the  extrapolated  field  representation.  An  alternate 
field  extrapolation  algorithm  will  be  stated,  which,  when  averaged  with  the  results  of  the 
first  algorithm,  cancels  the  effects  of  the  obliquity  factors  to  first  order.  The  two  estimates 
will  be  compared  to  monopole  and  dipole  secondary  source  contributions  from  Huygens’ 
principle.  It  will  then  be  shown  that,  for  range-dependent  media  (both  in  the  adiabatic 
approximation  and  for  coupled  mode  environments),  acoustic  calibration  takes  the  range- 
dependence  into  account  and  can  provide  accurate  extrapolations.  The  idea  of  reciprocal 
calibration  will  be  briefly  mentioned. 

Chapter  3  will  illustrate  some  of  the  mathematical  concepts  from  Ch.  2  through 
numerical  simulation.  Specifically,  it  will  explore  the  sensitivity  of  the  extrapolation 
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algorithm  to  one  fundamental  property  of  an  underwater  acoustic  environment  -  the 
“steepness”  of  energy  propagating  through  it.  It  will  be  shown  how  environments  that 
support  energy  with  high  grazing  angles  are  more  severely  impacted  by  the  higher  order 
terms  of  the  obliquity  factors.  Several  methods  of  comparing  extrapolated  and  “ground 
truth”  fields  will  be  presented.  Two  qualitative  comparisons  are  transmission  loss  plots  and 
phase  space  representations  (Argand  diagrams).  A  quantitative  measure  called  the 
matching  coefficient  will  be  introduced.  A  method  for  approximating  the  degradation  in 
the  matching  coefficient  based  on  approximations  of  the  obliquity  factors  will  be  detailed. 

Chapter  4  will  explore  the  sensitivity  of  the  extrapolation  algorithm  to  mismatch  in 
the  modeled  fields.  Mismatch  occurs  when  insufficient  and/or  incorrect  knowledge  of  the 
parameters  of  the  medium  that  govern  acoustic  propagation  are  used  as  input  to  the 
modeling  algorithm.  A  perturbational  approach  will  be  used  wherein  errors  in  the 
environmental  parameters  will  manifest  themselves  in  perturbations  to  the  normal  modes 
(eigenfunctions  and  eigenvalues),  and  the  previously  developed  normal  mode 
representations  of  the  fields  will  be  used  to  analyze  the  errors.  Mismatch  due  to  errors  in 
the  water  column  sound  speed  profile  and  some  sediment  parameters  (sound  speed, 
attenuation,  and  density)  will  be  considered.  Methods  for  approximating  the  degradation 
in  the  matching  coefficient  based  on  approximations  of  the  horizontal  wavenumber 
perturbations  (as  for  the  obliquity  factors)  will  be  detailed. 

Chapter  5  will  begin  by  exploring  the  sensitivity  of  the  extrapolation  algorithm  to 
the  spacing  of  transmitters/receivers  along  the  reference  arrays.  It  will  be  shown  how  a 
Fourier  theory  approach  can  be  used  to  analyze  how  increasing  the  spacing  can  degrade 
the  performance  of  the  algorithm  when  trapezoidal,  i.e.,  simple  summation,  numerical 
integration  (quadrature)  techniques  are  used.  It  will  be  shown  how  the  use  of  an  inner 
product  quadrature  formula  (I.P.Q.F.),  which  explicitly  utilizes  the  mode  functions  of  the 
modeled  field,  can  be  used  to  estimate  the  integral.  This  technique  requires  that  the 
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number  of  transmitters/receivers  along  the  reference  arrays  be  greater  than  or  equal  to  the 
number  of  modes  that  contribute  significantly  to  the  fields  for  that  environmental  scenario. 
It  will  also  be  shown  how  the  I.P.Q.F.  can  be  modified  in  order  to  cancel  the  obliquity 
factors. 

Chapter  6  will  then  summarize  the  dissertation  and  provide  suggestions  for  future 


work. 


Chapter  2:  Acoustic  Field  Extrapolation  -  Theory 


2.1.  Introduction  to  Chapter  2 

This  chapter  will  outline  the  basic  theory  behind  the  procedure  for  extrapolating 
acoustic  fields  in  the  ocean  through  the  combined  use  of  measured  acoustic  data  and 
acoustic  field  modeling  methods.  An  idealized  geometry  where  the  propagation  conditions 
are  known  exactly  will  be  described.  The  idea  of  “acoustic  calibration”  will  then  be 
discussed,  as  will  Huygens’  principle  and  its  relevance  to  the  field  extrapolation  algorithm. 
Mathematical  formulations  for  the  ideal  scenario  will  be  stated,  and  an  expression 
involving  a  combination  of  field  estimates  with  different  “obliquity  factors,”  which  can  be 
seen  to  approximate  contributions  from  monopole  and  dipole  secondary  sources,  will  be 
presented.  Range-dependent  scenarios  will  then  be  discussed. 

2.2.  Problem  Discussion 

A  situation  where  acoustic  field  prediction  is  required  is  shown  in  Figure  2.1.  We 
want  to  predict  the  field  at  the  arbitrarily  placed  hypothetical  receiver  due  to  the  arbitrarily 
placed  hypothetical  source.  It  is  known  that  for  short  ranges  in  the  vicinity  of  the 


Transmit  Receive 

Array  Array 


Figure  2.1.  Scenario  of  extrapolation  problem. 
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hypothetical  source  and  receiver  that  reliable  knowledge  about  the  environmental 
conditions  that  affect  acoustic  propagation  can  be  obtained,  and,  therefore,  that  reliable 
estimates  of  the  acoustic  fields  can  be  calculated  through  numerical  modeling.  It  is  also 
known,  however,  that  a  relatively  large  region  between  the  source  and  receiver  has 
unknown  or  difficult  to  measure  acoustic  propagation  properties,  leading  to  difficulties  in 
predicting  acoustic  fields  across  this  region  through  numerical  modeling  methods  alone. 

Now  we  place  vertical  arrays  of  sources  and  receivers  near  the  boundaries  of  the 
“unknown”  region  so  that  by  transmitting  and  receiving  acoustic  energy  through  this  part 
of  the  medium  it  is  “calibrated.”  This  situation  is  much  like  that  seen  in  the  HAP  algorithm 
in  Section  1.3,  and  that  development  should  give  us  a  clue  as  to  how  to  proceed  in 
predicting  the  acoustic  fields  between  hypothetical  sources  and  receivers  that  are  located 
outside  the  calibrated  region.  In  the  next  section,  we  will  examine  Huygens’  principle 
which  will  give  us  further  guidance  in  how  to  proceed. 

2.3.  Huygens’  Principle 

The  Dutch  physicist  Christian  Huygens  (1629-1695)  presented  his  theory  on  the 
nature  of  light  in  1678.  Huygens  postulated  that  if  the  position  of  a  wavefront  was  known 
at  time  1  =  0,  then  any  point  on  the  wavefront  could  be  considered  as  the  location  of  a 
“secondary  source”  which  could  produce  secondary  spherical  wavelets.  At  some  later  time 
t  =  tj ,  the  position  of  the  wavefront  would  be  the  tangent  surface  to  these  secondary 
wavelets.  This  idea  is  illustrated  in  Fig.  2.2.  Huygens’  principle  can  be  applied  to  all  wave 
phenomena  [Halliday  and  Resnick],  including  acoustic  waves. 

An  interesting  application  of  Huygens’  principle  from  the  field  of  optics  serves  as 
a  motivation  for  our  underwater  acoustic  application.  Imagine  the  situation  pictured  in 
Fig.  2.3  with  an  infinite  plane  screen  and  a  finite  aperture.  A  point  source  at  P2  illuminates 
the  aperture  with  the  field 


Figure  2.3.  Infinite  plane  screen  with  finite  aperture. 
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U(P  )  =  - ,  (2.1) 

'*21 

where  A  is  an  arbitrary  amplitude  factor,  is  the  length  of  the  vector  pointing  from  P2 
to  Pj ,  and  k  =  2k/?i  is  the  wave  number  with  X  the  wavelength  of  the  assumed 
monochromatic  light.  The  so-called  Huygens-Fresnel  principle  (an  interpretation  of  the 
Fresnel-Kirchhoff  diffraction  formula)  [Goodman]  states  that  the  field  at  Pq  can  be 
approximated  by 

v(p„)sffmp,f—ds, 

V  ^01 

^  (2.2) 

where  S  denotes  the  two-dimensional  area  of  the  aperture  and 


7^ ''21 


cos  (n,  Ff,,)  -  cos  (n,  f,,) 


=  A^U(Pj)0(n,  Fqj,  r2i). 


where  N  =  (jX)  cos  (a,  B)  denotes  the  cosine  of  the  angle  between  the  two  vectors, 
and  9  will  be  referred  to  as  an  “obliquity  factor.”  The  rest  of  the  integrand  in  Eq.  2.2  (apart 
from  U'(Pj))  can  be  thought  of  as  the  free-space  Green’s  function  for  the  region  on  the 
right  of  the  screen. 

Eq.  2.2  implies  that  “the  field  at  Pq  arises  from  an  infinity  of  fictitious  ‘secondary’ 
point  sources  located  within  the  aperture  itself.”  As  mentioned  above,  0  is  called  an 
obliquity  factor  and  applies  a  scaling  to  each  secondary  source  that  is  dependent  on  the 
directions  of  the  true  point  source  and  the  observation  point  relative  to  the  position  of  the 
given  secondary  source.  This  can  be  thought  of  as  applying  a  “directivity  pattern”  to  each 
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secondary  source.  We  will  draw  analogies  to  this  situation  in  our  acoustic  field 
extrapolation  derivation. 

2.4.  An  Idealized  Geometry  and  Acoustic  Calibration 

Figure  2.4  is  a  schematic  of  the  idealized  geometry  that  will  be  used  for  the  initial 
development  of  the  theory  behind  the  acoustic  calibration  algorithm.  Ultimately,  we  will 
be  interested  in  the  acoustic  pressure  p  =  p{r^,  at  the  receiver  due  to  a  source  emitting 
a  cw  signal  of  frequency  f  from  the  coordinates  (0,  Zq)  .  Our  ideal  ocean  will  be  assumed 
to  be  a  time-invariant  horizontally  stratified  acoustic  waveguide  of  depth  h  meters 
characterized  by  a  range-independent  sound  speed  profile  c(z).  The  ocean  surface  is 
modeled  as  a  planar  free  (or  pressure  release)  surface,  and  the  ocean  bottom  is  modeled  as 
a  perfectly  rigid  planar  surface  (see  Sec.  1.2.1  for  a  discussion  of  boundary  conditions). 

For  this  development,  the  water  column  will  be  assumed  to  be  spanned  by  an  ideal 
transmitting  array  at  range  and  an  ideal  receiving  array  at  range  as  shown  in  Fig.  2.4. 
These  reference  arrays  are  ideal  in  the  sense  that  sound  can  be  transmitted  from  any  point 
along  the  transmit  array  such  that  0<z^<h,  and  sound  can  be  received  at  any  point 
along  the  receive  array  such  that  0<Zy<h,  i.e.,  the  arrays  are  not  limited  by  having  fixed, 
discrete  transducer  locations.  This  allows  for  continuous  integration  of  functions  along  the 
arrays.  For  the  more  realistic  situation  of  fixed  transducer  locations,  numerical  integration 


Figure  2.4.  Idealized  geometry  for  acoustic  calibration  algorithm. 
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(quadrature)  techniques  can  be  used.  Also,  the  assumption  will  be  made  that  the  source, 
receiver,  and  reference  arrays  are  all  in  a  common  vertical  plane. 

This  acoustic  field  extrapolation  algorithm  is  based  on  the  idea  of  acoustic 
calibration  of  the  medium  as  was  the  HAP  algorithm  in  Sec.  1.3.  A  given  region  of  ocean 
is  calibrated  by  transmitting  a  cw  signal  of  frequency  f  from  all  depths  Q<z^<h  along 
the  reference  transmit  array  and  recording  the  resulting  complex  pressure  field  as  a 
function  of  depth  along  the  reference  receive  array  for  all  depths  0  <  <  /i .  This  process 

essentially  gives  us  a  spatial  impulse  response,  or  Green’s  function,  describing  the 
propagation  of  sound  from  to  r^,  G(r^,  z^;r^,  z^) .  It  will  be  assumed  that  the  reference 
transmit  array  and  the  hypothetical  source  operate  at  the  same  frequency  f,  which  will 
then  be  implicit  for  all  expressions  of  pressure  and  Green’s  functions  to  follow. 

The  idea  now  is  to  use  this  “calibration”  data  to  predict  what  the  acoustic  pressure 
would  be  at  the  receiver  at  coordinates  (r^,  z^)  due  to  the  source  at  coordinates  (0,  Zq)  . 
Once  the  calibration  (Green’s  function)  data  have  been  collected,  the  algorithm  proceeds 
as  follows:  1)  estimate  a  value  related  to  the  pressure  field  (to  be  described  later)  at  the 
transmit  reference  array  (r  =  r^)  due  to  the  source  using  local  environmental  information 
at  the  array  and  an  acoustic  field  modeling  method,  2)  combine  the  results  of  step  1)  with 
the  measured  Green’s  function  via  a  specialization  of  Huygens’  principle  to  obtain  an 
estimate  of  the  pressure  field  along  the  receive  reference  array  (r  =  r^),  3)  estimate  a 
value  related  to  the  Green’s  function  (to  be  described  later)  between  the  receive  reference 
array  and  the  receiver  (between  and  r^)  using  local  environmental  information  at  the 
array  and  an  acoustic  field  modeling  method,  and  4)  combine  the  result  of  step  2)  with  the 
result  of  step  3)  via  a  specialization  of  Huygens’  principle  to  obtain  an  estimate  of  the 
pressure  at  the  receiver. 

The  calibration  range  is  assumed  to  be  much  greater  than  either  of 

the  “extrapolation”  ranges,  or  r^2  =  ^c~^b-  behind  the 
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algorithm  is,  as  mentioned  in  Sec.  2.2,  that  the  calibration  procedure  will  measure  the 
acoustic  transmission  properties  of  the  medium  between  the  arrays  and  incorporate  this 
knowledge  into  the  field  estimate,  removing  the  effects  of  medium  variability  that  would 
be  difficult  to  determine  using  only  environmental  measurements. 

2.5.  Mathematical  Statement  of  Acoustic  Field  Extrapolation  Algorithm 

For  this  part  of  the  development  it  is  assumed  that  the  propagation  properties  of  the 
medium,  i.e.,  the  range-independent  sound  speed  profile  c(z)  and  subsequently  the  normal 
mode  decomposition,  are  known  perfectly  for  the  modeling  steps.  If  this  were  the  case  in 
true  applications,  there  would  be  no  advantage  to  using  the  extrapolation  algorithm  since 
one  of  the  other  propagation  simulation  methods  mentioned  in  Ch.  1  could  be  used  to 
predict  the  acoustic  field  over  the  entire  range.  This  assumption,  however,  is  helpful  in 
developing  the  theory.  The  more  realistic  cases  of  range-dependent  propagation  and 
mismatch  between  the  assumed  properties  and  the  true  properties  will  be  dealt  with  in 
later  sections. 

A  normal  mode  representation  will  be  used  for  both  the  measured  (calibration)  and 
modeled  fields.  It  should  be  noted  that  this  representation  is  not  a  requirement  of  the 
theory  for  this  part  of  the  development,  i.e.,  the  modeled  fields  could  be  generated  using 
any  of  a  number  of  acoustic  field  modeling  methods.  This  representation  is  useful, 
however,  in  showing  the  mathematical  relationships  between  the  actual  and  extrapolated 
fields. 

2.5.1.  Field  Extrapolation  to  the  Reference  Receive  Array 

Assuming  that  the  medium  has  been  acoustically  calibrated  between  and 
(giving  us  G{r^^,  ),  the  first  step  in  predicting  the  field  at  (r^,  z^)  is  to  extrapolate 

a  modeled  field  (due  to  the  source  at  (0,  Zq)  )  along  the  reference  transmit  array  at  range 
to  the  receive  array  at  range  As  mentioned  above,  this  can  be  accomplished  by 
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combining  the  modeled  field  with  the  Green’s  function  data  via  a  specialization  of 
Huygens’  principle.  One  way  to  do  this  is  to  define  [Rouseff  (Feb.  1993)] 


Zf)  =  zjdz^  =  , 


(2.4) 


where  the  tilde  denotes  that  this  is  an  extrapolated  value,  the  prime  corresponds  to  this 
particular  method  of  combining  the  information,  Nj  is  a  scaling  constant  to  be 
determined,  and  p(f^,  z^)  is  a  modeled  field  along  the  reference  transmit  array  due  to  the 
hypothetical  source.  This  is  a  modification  of  Eq.  2.2  where  we  are  now  integrating  in  one 
dimension,  the  measured  Green’s  function  G(r^,  zj  takes  the  place  of  the  free  space 
Green’s  function  and  pir^jZ^  takes  the  place  of  U'(Fi)-  Note  that  the 

“incident”  field,  p{r^,  is  not  modified  as  its  counterpart  is  in  Eqs.  2.2-3. 

Using  the  range-independent  normal  mode-based  expression  for  underwater 
acoustic  fields  (Eq.  1.10),  we  can  express  the  modeled  field  p{r^,  z^  as 


p{ra,  z^)  =  p{r^,  zj  = 


87ir„ 


1/2 


-1/2  ‘L’-a 


p(Zo) 


(2.5) 


The  caret  indicates  that  this  is  a  modeled  pressure  field  based  on  local  environmental  data. 
Since  the  assumption  of  perfect  knowledge  about  the  environment  has  been  made,  the 
estimate  is  equal  to  the  true  pressure  field  (the  approximate  equality  is  used  due  to  the 
Hankel  function  approximation  from  Eq.  1.9).  Assuming  that  r^,  the  Green’s  function 
can  also  be  represented  using  a  normal  mode  decomposition  by 


G(r^,  zJ  = 


^  ■  (2.6) 


Combining  Eqs.  2.4-6  and  taking  the  integral  inside  the  summations  we  find  that 


j  ,  _ _ i _ ]_  ^ 


,5.  P  ,-1/2!  j 

J-— ^?r; — ■ 


Since  the  integral  at  the  end  of  Eq.  2.7  is  a  Kronecker  delta  function  5^^  (due  to  the 
orthonormality  property  of  the  mode  functions  -  see  Sec.  1.2.1),  the  double  summation 


collapses  and 


The  true  pressure  field  at  the  receive  reference  array  would  be 


-1/2  '^''6 


'’‘''O' -  |_8^J  e 


Note  that  the  differences  between  Eqs.  2.8  and  2.9  occur  in  the  range  scaling  in  front  of  the 
summation,  and  the  power  of  the  horizontal  wavenumber  inside  the  summation  (-1  vs. 
-1  /2 ).  N^'  can  now  be  chosen  so  as  to  make  z^)  a  better  approximation  to  pir^^,  z^) 
[Rouseff  (Feb.  1993)].  By  setting 


^  ^  \Sn^refa(’-b-0 

"  I  ru 


(2.10) 


where  is  a  reference  wavenumber  (see  Eq.  1.13  and  accompanying  discussion), 
Eq.  2.4  becomes 
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—1/2 

Note  that  the  extra  has  been  combined  with  the  in  order  to  form  a  modal- 

dependent  multiplier.  This  multiplier  will  be  discussed  in  the  next  section. 

Care  must  be  taken  in  applying  the  normalization  constant  Nj  in  practice, 
especially  if  absolute  field  level  is  the  desired  result.  This  is  due  to  the  fact  that  some 
numerical  modeling  codes  use  different  conventions  for  normalizing  their  output  fields 
(see  Eq.  1.10). 


2.5.2.  The  Obliquity  Factor 

The  expressions  for  the  “ground  truth”  field  p(r^,  z^)  and  the  extrapolated  field 
p'(r^,  Zj)  due  to  the  hypothetical  source  differ  only  in  the  modal-dependent  multiplier 


(2.12) 


where  we  will  call  0^'  an  “obliquity  factor”  [Rouseff  (April  1993)].  In  referring  back  to 
the  Huygens-Fresnel  principle  (Eqs.  2.2-3),  which  was  the  motivation  for  our  acoustic 
field  extrapolation  formulation,  we  see  that  the  obliquity  factor  (which  depended  on  the 
source-aperture -receiver  geometry)  was  included  in  the  integration  as  a  modifier  to  the 
incident  field.  Here,  however,  the  incident  (modeled)  field  was  used  unaltered,  and  the 
obliquity  factor  appears  in  our  estimate  of  the  observed  field.  Indeed,  if  a  normal  mode 
formulation  is  used  for  the  modeled  field,  a  modified  scaling  of  each  mode  could  be 
applied  to  account  for  the  obliquity  factor.  For  now,  however,  we  will  assume  a  general 
formulation  in  that  no  modifications  of  this  kind  will  be  made  to  the  incident  field.  This 
allows  for  the  use  of  any  type  of  acoustic  field  modeling  algorithm  in  generating  the 
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modeled  fields. 

In  addition,  the  obliquity  factor  takes  on  a  geometric  interpretation  as  it  does  in  the 
optics  application.  Keeping  in  mind  that  the  approximate  grazing  angle  for  each  mode  can 
be  given  by 


Xm  ~ 


we  can  say  that  for  each  mode  in  the  estimated  field 


(2.13) 


(cosx^)-i/2.  (2.14) 

Since  lower  order  modes  propagate  with  grazing  angles  close  to  horizontal,  we  can  say 
that  6^'  =  1  for  these  modes.  Progressively  steeper  grazing  angles,  however,  lead  to  an 
overweighting  of  the  higher-order  modes  compared  to  the  true  pressure  field.  For 
environments  where  high-order  modes  propagate  over  long  distances,  this  can 
significantly  degrade  the  estimate  (examples  will  be  shown  in  Ch.  3).  Figure  2.5  shows  a 
plot  of  the  obliquity  factor  vs.  grazing  angle.  The  gray  line  shows  6^'  =  1  for  reference. 


Figure  2.5.  Obliquity  factor  vs.  grazing  angle. 
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The  obliquity  factor  tends  to  infinity  as  the  grazing  angle  tends  to  90' 


2.5.3.  An  Alternate  Extrapolation  Formula  and  Obliquity  Factor 

We  can  form  an  alternate  extrapolation  formula  by  defining  [Rouseff  (April  1993)] 


h  A  ✓  X-, 


(2.15) 


Here  the  secondary  sources  take  on  the  amplitude  and  phase  of  the  range  derivative  of  the 
incident  field.  In  practice  the  derivative  of  the  modeled  quantity  can  be  obtained  through 
discrete  difference  methods  or  direct  computation.  Differentiating  Eq.  2.5  gives 


1  -1/2  ’'a 


~  r  * 

~  Lsti. 


(2.16) 


We  will  neglect  the  range  term  to  the  -3/2  power  for  this  analysis  since  it  is  dominated  by 
the  other  range  term  as  range  increases.  So  doing  will  greatly  simplify  the  following 
results. 

Inserting  this  result  into  Eq.  2.15  and  taking  advantage  of  the  orthonormality  of  the 
mode  functions  as  before  gives 


(2.17) 


As  was  done  for  N^'  and  p'ir^,  Zj),  NJ'  can  be  chosen  so  as  to  make  p"(r^,  z^)  a 
better  approximation  to  p{r,,  Zi)  by  setting 


k  r 

L  I  J 


(2.18) 
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Grazing  Angle  (deg) 

Figure  2.6.  Alternate  obliquity  factor  V5.  grazing  angle. 

which  leads  to 


P"(r„z,)  = 


S>Ttr, 


n  1/2 


ref. 


(2.19) 


In  this  case  the  obliquity  factor  is 


0  "  = 

m 


1/2 


~  (cosx^) 


(2.20) 


which  leads  to  an  underweighting  of  the  higher-order  modes  compared  to  the  true  pressure 
field.  Figure  2.6  is  a  plot  of  this  alternate  obliquity  factor  vs.  grazing  angle.  This  obliquity 
factor  tends  to  zero  as  the  grazing  angle  tends  to  90° . 

The  Huygens-Fresnel  principle  (an  interpretation  of  the  Fresnel-Kirchhoff 
diffraction  formula)  is  an  approximate  solution  to  the  plane  screen  diffraction  problem  that 
is  derived  from  the  integral  theorem  of  Helmholtz  and  Kirchhoff  (see  [Goodman]  and 
[Porter  et  al.  (1992)]).  This  theorem  can  be  stated  as 
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5  (2.21) 

where  5  is  a  closed  surface  surrounding  the  point  Pq,  G  is  the  Green’s  function,  and  the 
derivatives  are  normal  to  the  surface.  An  interpretation  of  this  formula  from  the  standpoint 
of  Huygens’  principle  is  that  the  field  at  Pq  is  the  result  of  both  monopole  secondary 
sources  that  radiate  the  Green’s  function  and  dipole  secondary  sources  that  radiate  the 
normal  derivative  of  the  Green’s  function  [Porter  (1994)]. 

The  extrapolation  algorithm  of  Eq.  2.15  can  be  thought  of  as  an  approximation  to 
the  monopole  contribution  from  the  Helmholtz-Kirchhoff  integral,  since  we  can  think  of 
the  secondary  sources  having  the  amplitude  and  phase  of  the  spatial  derivative  of  the 
incident  field,  and  radiating  the  Green’s  function.  The  extrapolation  algorithm  of  Eq.  2.4 
can  be  thought  of  as  an  approximation  to  the  dipole  contribution  from  the  Helmholtz- 
Kirchhoff  integral  in  the  sense  that  the  secondary  sources  take  on  the  amplitude  and  phase 
of  the  incident  field  itself  rather  than  its  spatial  derivative.  In  general,  the  spatial 
derivatives  of  acoustic  fields  are  difficult  to  measure,  and  since  for  the  extrapolation 
algorithm  the  Green’s  function  is  measured  data,  we  consider  the  approximation 

(2-22) 

«  m 

where  (p  is  a  geometry  dependent  multiplier.  It’s  form  is  similar  to  the  multiplier  seen  in 
Eq.  2.16,  and  is  encompassed  in  the  obliquity  factor  0^'. 

2.5.4.  The  Averaged  Obliquity  Factor 

We  now  have  two  ways  to  synthesize  our  estimated  field,  one  that  overweights  and 
one  that  underweights  contributions  according  to  grazing  angle.  While  the  Helmholtz- 
Kirchhoff  integral  implies  that  a  differencing  operation  should  take  place  between  the 
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Figure  2.7.  Averaged  obliquity  factor  V5.  grazing  angle. 

monopole  and  dipole  contributions  in  order  to  determine  the  observed  field,  we  have 
artificially  scaled  the  two  contributions  so  that  they  are  each  better  approximations  to  the 
observed  field.  In  this  case,  averaging  the  two  components  is  more  appropriate,  i.e.. 


2 


(2.23) 


Since  Eqs.  2.11  and  2.19  differ  only  in  the  obliquity  factors,  we  can  think  of  their  average 
as  having  a  new  mode-dependent  obliquity  factor  6^.  Fig.  2.7  shows  a  comparison  of  the 
two  previously  plotted  obliquity  factors  (dashed  lines)  with  the  averaged  one  (solid  line) 
vs.  grazing  angle. 

As  an  example,  consider  a  mode  with  a  grazing  angle  of  60  degrees.  Its 
contribution  to  the  estimate  p'(r^,  z^)  would  be  -1.41  times  what  its  contribution  to  the 
true  pressure  field  p{r^^,  zf)  would  be,  while  its  contribution  to  the  estimate  p'\r^^,  zf) 
would  only  be  ~0.71  times  what  it  should  be.  Averaging  the  two  estimates  results  in  a 
contribution  of  -1.06  times  its  proper  contribution.  The  averaged  obliquity  factor  is  closer 
to  1  than  either  of  the  other  two  for  all  but  the  steepest  grazing  angles,  but  it  still  tends  to 
infinity  as  the  grazing  angle  moves  toward  90  degrees. 
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One  way  to  analyze  the  averaging  effect  is  to  consider  the  as  perturbations 
from  the  reference  wavenumber  [Rouseff  (April  1993),  Porter  et  al.  (1993),  Fox  et 
all 

-ej.  (2-24) 

where  is  assumed  small.  Expanding  the  obliquity  factors  in  a  Taylor  series  about 
=  0  gives 


K  =  <1 

Averaging  the  right  hand  sides  results  in  a  cancellation  of  the  first  order  terms.  The  field 
estimate  p  (r^,  z^)  is,  therefore,  said  to  be  “invariant”  to  first  order  errors  introduced  by 
the  obliquity  factors  [Porter  et  al  (1993)]. 

It  may  be  surmised  from  Fig.  2.7  that  the  higher  order  terms  from  Eq.  2.23  become 
significant  for  the  steeper  grazing  angles.  If  we  continue  the  Taylor  series  expansion  out  to 
second  and  third  order  terms,  we  obtain 


(2.26) 


Since  >  0 ,  averaging  the  right-hand  sides  of  Eq.  2.26  shows  that  second  and  third  order 
terms  are  not  canceled  and  are  positive.  This  accounts  for  the  upturn  in  the  averaged  curve 
of  Fig.  2.7.  Note  that  the  Helmholtz-Kirchhoff  integral  is  an  exact  formulation.  Our  errors 
are  due  to  the  facts  that  we  are  integrating  in  one  dimension  along  the  linear  reference 
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array  instead  of  over  a  closed  surface  and  that  we  have  used  the  measured  Green’s 
function  instead  of  its  spatial  derivative  in  the  dipole  contribution. 

2.5.5.  Field  Extrapolation  to  the  Receiver 

The  value  z^)  can  now  be  extrapolated  further  in  order  to  predict  the  pressure 
field  at  the  receiver  located  at  coordinates  (r^,  z^)  [Rouseff  (Feb.  1993)].  In  extrapolating 
out  to  the  reference  receive  array,  environmental  data  local  to  the  reference  transmit  array 
were  used  as  input  to  an  acoustic  field  modeling  routine  to  predict  a  pressure  field  and  its 
range  derivative  at  the  reference  transmit  array  (due  to  the  hypothetical  source  at 
coordinates  (0,  Zq)  ).  These  estimates  were  then  combined  with  the  measured  Green’s 
function  data  in  order  to  calculate  the  extrapolated  field.  Extrapolating  from  to  is  a 
dual  problem  of  the  first  case.  Now  the  Green’s  function  between  r  =  and  r  =  r^, 
G(r^,  Zf) ,  and  its  range  derivative  at  r  =  must  be  modeled  based  on 
environmental  data  at  the  reference  receive  array,  and  then  combined  with  p{r^,  z^ . 

Specifically,  we  will  now  define 

h 

~p\rc,  z^)  =  N,;yp{r^,  z^)G{r^,  z^;r^,  Zf^dz^  =  .  (2.27) 

0 

Using  the  approximation  p(/-^,  z^)  =  p(r^,  z^),  the  math  follows  much  as  before,  and  we 
can  set 


’Snkrefb(r,-r,,) 


-11/2 


(2.28) 


We  also  define 
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p\rc,z^)  =  Nj;'\~p{rb,  z^) 


h')' 

dr. 


dzi 


(2.29) 


and  set 


v  = 


Sni- 


k  fV 
ref  c  J 


n  1/2 


(2.30) 


The  obliquity  factors  work  out  as  before,  so  our  final  estimate  of  the  pressure  at  the 
receiver  is 


P(''c’  ^c) 


~p\rc,  z^)  +  p"{r^,  z,) 
2 


(2.31) 


The  flow  diagram  in  Fig.  2.8  gives  an  overview  of  the  extrapolation  algorithm. 
Note  that  the  part  of  the  algorithm  that  deals  with  extrapolating  from  the  receive  reference 
array  to  the  hypothetical  receiver  could  also  be  used  alone  in  the  case  where  there  was  an 
actual  source,  there  was  no  transmit  reference  array,  and  the  field  generated  by  the  source 
was  measured  along  the  receive  reference  array.  In  other  words,  measuring  the  field  along 
the  array  would  allow  one  to  extrapolate  that  field  for  hypothetical  receiver  locations 
beyond  the  array. 


2.5.6.  Adiabatic  Propagation  and  Acoustic  Field  Extrapolation 

We  now  consider  the  situation  where  the  medium  varies  slightly  with  range  and  the 
adiabatic  approximation  adequately  describes  the  pressure  fields.  This  section  will 
consider  only  the  part  of  the  extrapolation  algorithm  that  estimates  the  field  along  the 
receive  reference  array  due  to  a  hypothetical  source  outside  the  calibration  region. 

First  we  will  assume  that  the  medium  is  adiabatic  between  the  reference  arrays  and 
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Figure  2.8.  Flow  diagram  of  extrapolation  algorithm. 


range-independent  outside  this  region  (much  as  is  pictured  in  Fig.  2.1).  We  can,  therefore, 
write  the  true  pressure  along  the  receive  reference  array  as 


-  i  1 

_87trJ  p(zq) 


X 


^ 


(2.32) 


keeping  in  nrind  that  vi/^(r,  z)  =  z)  and  ^^(r)  =  ^J.r^)  for  0  <  r  <  due  to  the 

range-independence  in  the  region  near  the  hypothetical  source.  We  can  also  write  the 
measured  Green’s  function  as 
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G{rb,z^-,r^,  z„)  = 


-il/2 


-rj 


P(^«) 


^  -1/2 

^ 


(2.33) 


Applying  the  specialization  of  Huygens’  principle  as  for  the  previous  case  (using 
Eqs.  2.4-5,  2.10,  and  2.33),  we  find  that 


P'(fb’  h')  = 


Sttr 


bJ 


1/2  I 


Pizj 


X 


—1  /2 

XVm(''0’  ^ 


^ref 

1/2 


(2.34) 


The  only  difference  between  the  extrapolated  field  (Eq.  2.34)  and  the  true  field  (Eq.  2.32) 
is  the  modal-dependent  obliquity  factor,  as  was  the  case  for  the  ideal  range-independent 
case.  The  acoustic  calibration  step,  therefore,  takes  into  account  the  range-dependence  of 
the  medium  that  would  be  difficult  to  measure  by  direct  sampling  of  the  environmental 
properties.  The  same  technique  of  using  the  range  derivative  of  the  modeled  pressure  field 
at  the  reference  array  to  form  another  estimate  and  averaging  the  two  estimates  can  be 
used  to  ameliorate  the  effects  of  the  obliquity  factor. 

2.5.7.  Coupled  Modes  and  Acoustic  Field  Extrapolation 

We  now  will  show  how  the  acoustic  calibration  step  can  also  take  into  account 
more  severe  environmental  variability  to  aid  in  predicting  acoustic  fields.  For 
environments  where  the  adiabatic  approximation  does  not  hold,  we  can  use  the  coupled 
mode  formulation  for  underwater  acoustic  fields. 

In  order  to  keep  the  math  tractable,  we  will  work  with  the  problem  of  extrapolating 
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Figure  2.9.  Scenario  for  extrapolation  in  coupled  mode  environment. 


a  field  received  along  a  reference  receive  array  as  pictured  in  Fig.  2.9  (it  was  noted  earlier 
that  this  is  a  dual  problem  of  predicting  the  field  along  an  array  due  to  a  hypothetical 
source).  We  are  assuming  that  the  acoustic  energy  from  the  source  has  travelled  through  a 
part  of  the  medium  that  must  be  described  using  coupled  mode  theory,  and  that  the 
extrapolation  region  is  range-independent. 

Using  the  coupled  mode  theory  described  in  Sec.  1.2.3,  we  can  write  the  pressure 
field  measured  along  the  receive  array  as 

(2.35) 

m 

and  the  desired  pressure  field  at  the  hypothetical  receiver  as 

(2-36) 

m 

since  m  ~  ^  [Evans]  (see  Eq.  1.16).  We  write  the  estimated  Green’s  function  in  the 
region  near  the  reference  receive  array  as  (again  assuming  perfect  knowledge  of  the 
environmental  conditions  for  r>  r^) 
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G{r^,Z^',r^^,Zf)  =  G(r^,z^;ri^,Zi)  = 


i  1  X-  .  N  /  se-l/2 


which  is  our  usual  range-independent  estimate  with  perfect  knowledge  of  the 
environmental  parameters  that  affect  the  propagation. 

Combining  Eqs.  2.35  and  2.37  via  Eq.  2.27  and  taking  account  of  the 
orthonormality  of  the  modes  at  range  r  =  results  in 


,  r  i  .  rn  /  N  / 

'l>  iK(r  -r  ) 

t'c  J  m 


(2.38) 


Considering  the  asymptotic  expansion  of  Hl^^  Jr)  for  any  r  >  j , 


H\b,Jr) 


__2__-|’/2^'V« 


-1/2  .• 


(2.39) 


and  using  the  normalizing  constant  from  Eq.  2.28,  i.e.. 


^nKefb^Jc-Jb 


(2.40) 


we  find  that 


~PN^rc,Zc)  =  F - 

m  L  N,n 


KefV'^ 


(2.41) 


Eq.  2.41  is  the  now  familiar  result  of  the  true  pressure  field  with  modal-dependent 
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obliquity  factors  (compare  with  Eq.  2.36).  The  same  technique  of  averaging  estimates  can 
be  used  to  lessen  the  effect  of  the  obliquity  factor.  The  point  here  is  that,  as  for  the  case  of 
mild  range  dependence  where  the  adiabatic  approximation  could  be  used  to  describe  the 
calibrated  data,  the  calibration  step  can  also  effectively  take  into  account  environments  for 
which  coupled  mode  formulations  must  be  used  to  describe  the  propagation.  Once  this 
part  of  the  medium  has  been  calibrated,  reliable  field  extrapolations  can  take  place. 

We  also  briefly  note  that  known  range  dependencies  in  the  extrapolation  regions 
can  be  taken  into  account  by  the  modeled  fields  through  use  of  range-dependent  modeling 
schemes.  This  can  improve  the  fidelity  of  the  extrapolated  fields.  The  case  of  range- 
dependent  bathymetry  is  considered  in  [Rouseff  et  al.]. 

2.5.8.  Reciprocal  Calibration 

One  form  of  the  principle  of  acoustic  reciprocity  is  stated  simply  in  [Kinsler  et  al.} 
as  “if  in  an  unchanging  environment  the  locations  of  a  small  source  and  a  small  receiver 
are  interchanged,  the  received  signal  will  remain  the  same.”  An  informative  application 
which  utilizes  the  reciprocity  principle  can  be  found  in  [Dowling].  Applying  this  to  the 
situation  pictured  in  Fig.  2.4,  we  can  say  that  calibrating  by  transmitting  from  r  =  and 
receiving  at  r  =  is  the  same  as  calibrating  by  transmitting  from  r  =  and  receiving 
at  r  =  r^.  In  other  words,  fields  can  be  extrapolated  due  to  hypothetical  sources  either  in 
the  region  r<r^  or  r  >  using  the  one-way  calibration  data.  One  way  to  think  of  this  is 
to  write 


(2.42) 


where  T  represents  a  transpose  operation. 
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2.6.  Summary  of  Chapter  2 

This  chapter  has  outlined  the  basic  theory  behind  the  acoustic  field  extrapolation 
algorithm.  The  problem  of  predicting  acoustic  fields  over  long  ranges  that  include 
environments  with  unknown  or  difficult  to  determine  propagation  conditions  was 
discussed,  and  how  deploying  vertical  arrays  of  transmitters  and  receivers  may  aid  in  this 
effort.  Huygens’  principle  and  an  interpretation  of  the  Fresnel-Kirchhoff  diffraction 
formula  (the  Huygens-Fresnel  principle)  were  discussed  in  order  to  motivate  our 
extrapolation  formula.  An  ideal  range-independent  ocean  waveguide  was  described,  and 
the  idea  of  acoustic  calibration  using  reference  arrays  spanning  the  medium  in  depth  was 
introduced. 

Perfect  knowledge  of  the  propagation  conditions  was  assumed,  and  a  mathematical 
statement  for  this  extrapolation  was  presented  that  used  a  modification  of  the  Fresnel- 
Kirchhoff  diffraction  formula.  An  obliquity  factor  in  the  expression  for  the  extrapolated 
field,  analogous  to  that  seen  in  the  optics  example,  was  described.  By  forming  two 
versions  of  the  extrapolated  field  and  averaging  them  (one  based  on  an  estimate  of  the 
pressure  field  at  the  reference  transmit  array  and  the  other  on  its  range  derivative),  the  first- 
order  effects  of  the  obliquity  factor  were  negated.  The  two  estimates  were  likened  to 
monopole  and  dipole  secondary  source  contributions  from  Huygens’  principle.  It  was  then 
shown  how  the  field  can  be  extrapolated  to  a  receiver  in  the  vicinity  of  the  reference 
receive  array  using  similar  methods. 

Range-dependence  in  the  medium  (for  the  case  of  adiabatic  propagation  as  well  as 
full  mode  coupling)  was  discussed,  and  it  was  seen  that  the  extrapolation  formulas  stated 
for  the  range-independent  case  were  still  valid.  It  was  also  shown  how  the  roles  of  the 
reference  transmit  and  receive  arrays  can  be  interchanged  due  to  the  acoustic  reciprocity 
principle. 

The  following  chapter  will  illustrate  some  of  the  mathematical  relationships 
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expressed  in  this  chapter  through  numerical  modeling. 


Chapter  3:  Acoustic  Field  Extrapolation  -  Numerical  Modeling 

3.1.  Introduction  to  Chapter  3 

This  chapter  will  use  numerical  modeling  methods  to  illustrate  the  acoustic  field 
extrapolation  algorithm  discussed  in  Ch.  2.  These  examples  will  only  illustrate  the  part  of 
the  algorithm  that  extrapolates  the  acoustic  pressure  field  to  the  reference  receive  array 
(r  =  r^).  The  extrapolated  field  will  be  compared  to  a  “ground  truth”  estimate  of  the 
acoustic  field  (computed  using  perfect  knowledge  of  the  medium  and  a  conventional 
numerical  propagation  method)  along  the  array. 

The  normal  mode  computer  program  KRAKEN  [Porter  (1992)]  will  be  used  to 
synthesize  all  modeled  fields  in  these  examples.  The  first  examples  will  examine  the 
effects  of  the  obliquity  factor  on  the  extrapolation  algorithm  in  an  acoustic  waveguide  with 
perfectly  reflecting  boundaries  and  no  attenuation.  The  extrapolated  and  ground  truth 
fields  will  be  compared  qualitatively  through  the  use  of  transmission  loss  and  phase  space 
plots,  as  well  as  quantitatively  through  the  use  of  a  matching  coefficient.  Additionally,  a 
measure  that  characterizes  a  given  environment  called  the  obliquity  characterization 
number  will  be  defined,  and  its  relationship  to  the  matching  coefficient  will  be  examined. 
Finally,  the  effects  of  an  attenuating  sediment  layer  on  the  algorithm  will  be  illustrated. 

In  this  chapter  and  the  rest  of  the  dissertation,  examples  will  be  numbered  using 
the  convention  of  a  number  followed  by  a  letter,  where  the  number  is  the  chapter  number 
and  the  letter  (A,B,C,...)  differentiates  examples  within  that  chapter. 

3.2.  Obliquity  Factor  Effects:  Rigid  Bottom  (Example  3A) 

An  expression  was  derived  in  Sec.  2.5  for  the  extrapolated  pressure  field  at  the 
reference  receive  array,  p{r^,  z^,  due  to  a  source  at  position  (0,  Zq)  .  This  estimated  field 
is  the  result  of  averaging  two  extrapolated  fields  that  deviate  from  the  “ground  truth”  field 
by  modal-dependent  obliquity  factors.  By  averaging  these  two  estimates,  the  effect  of  the 
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obliquity  factors  is  cancelled  to  first  order.  This  example  will  utilize  a  simple 
environmental  scenario  to  illustrate  the  procedure. 

3.2.1.  Modal  Structure 

The  simulated  ocean  for  Ex.  3A  will  follow  the  idealized  geometry  of  Sec.  2.4.  We 
assume  a  pressure  release  ocean  surface  and  a  perfectly  rigid  bottom.  The  operating 
frequency  will  be  /  =  150  Hz.  The  depth  of  the  ocean  will  he  h  =  50  m ,  and  the  sound 
speed  profile  will  be  range-independent  and  also  constant  with  depth.  The  calibration 
range  will  be  =  20  km .  The  reference  transmit  and  receive  arrays  will  be  populated 
with  ideal  point  sources  with  spacing  Az^  =  Az^  =  1  m .  Figure  3. 1  is  a  schematic  for  this 
scenario. 

For  Ex.  3A  we  will  use  a  constant  sound  speed  c  (z)  =  1500  m/s.  KRAKEN 
calculates  10  propagating  modes  for  this  scenario,  and  the  mode  functions  are  plotted  in 
Fig.  3.2.  The  dashed  gray  lines  in  the  plots  represent  zero  pressure,  and  each  mode  is 
individually  normalized  by  the  maximum  amplitude  of  that  mode.  Also  note  how  each 
mode,  or  eigenfunction,  satisfies  the  boundary  conditions:  the  pressure  is  zero  at  z  =  0 
and  the  derivative  of  the  pressure  is  zero  at  z  =  h.  The  horizontal  wavenumber,  or 
eigenvalue,  of  each  mode  is  tabulated  (to  three  decimal  places)  along  with  the  approximate 
grazing  angle  (to  one  decimal  place)  in  Table  3.1.  Note  that  for  this  situation,  the  mode 

=  20  km - ► 


Figure  3.1.  Geometry  for  Exs.  3 A  and  3B. 


Depth  (m)  Depth  (m) 
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Figure  3.2.  Mode  functions  for  Ex.  3 A. 
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Table  3.1.  Eigenvalues  for  Ex.  3A. 


Mode  #  (m) 

Eigenvalue 

Grazing  Angle  (xj 
(deg.) 

1 

0.628 

2.9 

2 

0.621 

8.6 

3 

0.608 

14.5 

4 

0.589 

20.5 

5 

0.561 

26.7 

6 

0.525 

33.4 

7 

0.477 

40.5 

8 

0.416 

48.6 

9 

0.331 

58.2 

10 

0.196 

71.8 

functions  and  horizontal  wavenumbers  can  be  obtained  analytically  [Boyles]  and 
represented  as 


= 


7,2  (  {2m-l)K 

“I  h  2 


1/2 
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r2n  1/2 

.h. 


1/2  -| 


(3.1) 


for  m  =  1,2,  3, ... 

3.2.2.  Primary  Extrapolation 

To  simulate  calibration  of  the  medium,  we  numerically  model  the  fields  that  would 
be  incident  on  the  receive  reference  array  due  to  the  sources  along  the  transmit  reference 
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array,  giving  us  a  discrete  version  of  Gir^^,  Zfy\r^,  zj  ■  We  then  place  a  hypothetical  source 
=  3  km  from  the  reference  transmit  array  and  =  25  m  in  depth  and  numerically 
model  what  the  field  would  be  along  the  reference  transmit  array,  p(r^,  zj  ■  We  now  apply 
Eq.  2.4  in  order  to  extrapolate  the  field  to  the  reference  receive  array.  A  simple  trapezoidal 
rule  is  used  to  approximate  the  integral  (numerical  integration  issues  will  be  addressed  in 
Ch.  5). 

Figure  3.3  shows  comparisons  of  the  extrapolated  field  to  the  “ground  truth”  field 
which  is  calculated  by  numerically  modeling  the  field  along  the  receive  reference  array 
due  to  the  hypothetical  source.  Figure  3.3a  compares  the  transmission  loss  of  the  fields. 
Transmission  loss  is  defined  as 


TL(r,z)  =  -201og 


Pjr,  z) 


(3.2) 


where 


ikr 

0  e 

P  =  4i-r 


(3.3) 


is  the  expression  for  the  pressure  due  to  the  source  in  free  space.  Large  transmission  loss 
values  correspond  to  small  acoustic  pressures,  hence  the  inverted  nature  of  the  scaling  on 
the  y-axis. 

In  Fig.  3.3a  the  black  line  represents  the  ground  truth  transmission  loss  curve,  and 
the  gray  line  represents  the  extrapolated  estimate.  Note  how  the  extrapolated  transmission 
loss  curve  is  consistently  above  the  ground  truth  curve,  but  especially  so  near  the  peaks  in 
the  curve.  Equations  1.10-11  showed  how  the  expression  for  the  pressure  field  can  be 
thought  of  as  a  linear  sum  of  the  mode  functions  with  complex-valued  weights.  The 
extrapolated  field,  then,  can  be  thought  of  in  the  same  way,  except  with  the  weights 


Figure  3.3.  Ex.  3A: 
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multiplied  by  the  obliquity  factor.  Since  the  peaks  in  the  transmission  loss  curve  are  where 
the  contributions  to  the  pressure  field  tend  to  constructively  interfere  to  the  greatest 
degree,  the  overemphasis  caused  by  the  obliquity  factors  (9^')  is  most  evident  here.  Fig. 
3.3b  is  simply  the  difference  between  the  estimated  and  true  transmission  loss  curves  in 
Fig.  3.3a. 

No  information  about  the  phase  of  the  extrapolated  field  is  available  through 
examining  transmission  loss  only.  Figure  3.3c  is  a  so-called  “Argand  diagram”  [Uscinski]. 
It  shows  the  trajectory,  parameterized  by  depth,  of  the  complex-valued  pressure  fields  in 
phase  space.  The  outside  circle  represents  the  maximum  amplitude  contained  in  either 
field.  Again  the  ground  truth  is  the  black  line  and  the  extrapolated  is  gray.  Although  these 
diagrams  can  tend  to  be  quite  “busy,”  they  can  give  a  good  qualitative  feel  for  how  close  of 
a  coherent  “match”  the  two  fields  are. 

A  more  quantitative  measure  is  obtained  by  calculating  the  magnitude-squared  of 
the  “normalized  scalar  product”  [Goncharov  and  Voronovich]  of  the  two  pressure  fields. 
This  is  given  by 


\K\ 


2 


(3.4) 


2 

where  N  is  the  number  of  receivers.  For  this  work,  |^r|  will  be  referred  to  as  the 
“matching  coefficient.”  The  maximum  value  of  this  measure  is  unity,  and  is  attained  when 
the  two  fields  differ  only  by  a  constant  complex  multiplier.  For  the  situation  in  Fig.  3.3, 
\K\'^b  0.940. 
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3.2.3.  Alternate  and  Averaged  Extrapolations 

We  can  also  calculate  the  alternate  extrapolated  field  given  by  Eq.  2.15,  the  results 
of  which  are  shown  in  Fig.  3.4.  The  three  parts  of  Fig.  3.4  follow  the  same  pattern  as 
Fig.  3.3.  Note  that  now  the  amplitude  of  the  extrapolated  field  is  generally  underestimated 
in  comparison  with  the  ground  truth  field  due  to  the  attenuating  effect  of  the  obliquity 
factors  on  the  constituent  modal  amplitudes.  Here,  |.K’|  =0.936.  The  averaged 
extrapolated  field  (as  per  Eq.  2.23)  is  shown  in  Fig.  3.5.  The  partial  cancellation  of  the  first 
two  obliquity  factors  has  brought  the  matching  coefficient  up  to  |.^r|  =  0.982 . 

It  should  be  noted  that  as  the  location  of  the  hypothetical  source  is  varied  in  range 
and  depth,  the  nature  of  the  pressure  fields  produced  by  it  will  change.  For  example, 
placement  of  the  source  in  depth  determines  to  what  extent  each  mode  is  excited,  while 
placement  in  range  affects  the  phase  multiplier,  and,  hence,  the  constructive/destructive 
nature  of  the  mode  functions’  interaction  (see  Eq.  1.10).  Variation  of  the  source  position 
within  the  same  environmental  scenario,  then,  will  produce  small  fluctuations  in  the 
matching  coefficient  measure  for  two  reasons.  One  is  that  if  some  higher  order  modes  are 
not  excited  to  a  significant  extent,  this  will  lower  the  effects  that  their  obliquity  factors 
would  have  on  the  extrapolation.  The  second  is  that  the  matching  coefficient  is  somewhat 
sensitive  to  the  actual  shapes  of  the  pressure  waveforms.  Since  the  errors  in  the 
extrapolated  fields  vary  with  depth,  the  points  where  the  fields  happen  to  be  evaluated/ 
measured  (i.e.,  the  locations  of  receivers)  may  contribute  more  or  less  to  a  matching  error. 
This  contribution  to  the  fluctuations  will  decrease  as  the  number  of  receiver  points  is 
increased. 

Since  we  are  not  interested  in  these  incidental  fluctuations,  but  rather  the  overall 
performance  as  a  function  of  environment,  we  will  average  results  over  a  grid  of  source 
locations.  For  this  scenario,  the  hypothetical  source  location  was  varied  over  a  grid  of 
ranges  from  1  km  to  5  km  spaced  at  1  km  and  depths  from  5  m  to  45  m  spaced  at  5  m,  and 
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Alternate  field  extrapolation. 
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the  average  value  of  \K^  for  p{r,  z)  was  found  to  be  0.984  with  a  variance  of  7.545x10'^. 
The  range  and  depth  of  the  source  used  in  this  example  were  then  chosen  for  the  plots 
because  the  extrapolated  field  for  this  source  position  gave  a  matching  coefficient  close  to 
the  average  over  the  grid,  making  it  a  representative  example. 

3.3.  Obliquity  Factor  Effects:  Steeper  Grazing  Angles  (Example  3B) 

Example  3B  is  designed  to  increase  the  severity  of  the  effects  of  the  obliquity 
factors  in  comparison  with  Ex.  3A.  The  only  change  here  from  Ex.  3A  will  be  a  constant 
sound  speed  of  1570  m/s  instead  of  1500  m/s. 

3.3.1.  Modal  Structure 

KRAKEN  calculates  10  propagating  modes  for  Ex.  3B,  with  identical  mode 
functions  to  those  of  Ex.  3A  (see  Fig.  3.2).  This  may  be  surprising  at  first,  but  through 
examination  of  Eq.  3. 1  we  see  that  the  mode  functions  are  sections  of  sine  waves  that  must 
meet  certain  boundary  conditions  at  the  surface  and  bottom.  The  horizontal  wavenumbers 
(eigenvalues),  however,  have  changed  and  are  listed  in  Table  3.2.  Note  that  all  the  grazing 
angles  have  become  steeper,  leading  to  obliquity  factors  further  from  one  for  each  mode 
when  compared  with  those  of  Ex.  3A. 

3.3.2  Primary,  Alternate,  and  Averaged  Extrapolations 

Using  a  hypothetical  source  location  (r^,,  z^)  =  (2  km,  35  m)  ,  the  extrapolated 
vs.  ground  truth  field  comparisons  can  be  found  in  Figs.  3.6-8  with  analogous  ordering  to 
those  found  in  Ex.  3A.  The  matching  coefficients  are,  respectively,  0.848,  0.782,  and 
0.944.  For  the  same  grid  of  source  locations  as  in  Ex.  3A,  the  average  value  of  the 
matching  coefficient  was  0.945  with  a  variance  of  2.7 19x10''^.  Note  how  the  increase  in 
the  grazing  angles,  and  consequently  the  distance  of  the  obliquity  factors  from  unity,  has 
degraded  the  quality  of  the  extrapolations  in  comparison  with  Ex.  3A. 
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Table  3.2.  Eigenvalues  for  Ex.  3B. 


Mode  #  (m) 

Eigenvalue 

Grazing  Angle  (xj 
(deg.) 

1 

0.599 

3.0 

2 

0.593 

9.0 

3 

0.579 

15.2 

4 

0.559 

21.5 

5 

0.530 

28.1 

6 

0.491 

35.1 

7 

0.440 

42.9 

8 

0.372 

51.7 

9 

0.274 

62.8 

10 

0.064 

83.9 

3.4.  Obliquity  Characterization  Number 

At  this  point,  we  desire  some  simple  quantity  that  can  be  calculated  before  any 
fields  are  measured  or  modeled  in  order  to  aid  us  in  determining  to  what  extent  the 
acoustic  field  extrapolation  algorithm  will  be  affected  by  the  obliquity  factors  for  a  given 
underwater  environment.  Since  for  this  work  our  “goodness”  measure  is  the  matching 
coefficient,  we  will  use  that  as  our  starting  point. 

We  note  that  the  summations  over  the  discretely  sampled  points  in  depth  in  the 
definition  of  the  matching  coefficient  (Eq.  3.4)  can  be  thought  of  as  approximations  to  an 
integral  formulation,  i.e.. 
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where  the  integrals  are  over  the  length  of  the  array.  We  now  substitute  p(z)  (the 
extrapolated  field)  for  p  l(z)  and  p{z)  (the  ground  truth  field)  for  p2{z).  Noting  that  we  can 
write  the  fields  as 


Pte)  = 

m 

(3.6) 

where  the  6^  are  the  obliquity  factors,  we  substitute  into  the  second  part  of  Eq.  3.5  and 
obtain 
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(3.7) 


Interchanging  the  order  of  summation  and  integration  in  the  numerator  and 
expanding  the  magnitude  squared  operators  while  interchanging  the  order  of  summation 
and  integration  in  the  denominator  results  in 
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(3.8) 
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Assuming  real-valued  mode  functions  and  unity  density  profiles,  the  inner  products  of  the 
mode  functions  become  Kronecker  delta  functions  6^^  and  collapse  the  double 
summations.  Also  assuming  that  only  the  real-valued  parts  of  the  horizontal  wavenumbers 
are  used  to  calculate  the  obliquity  factors  and  expanding  the  magnitude  squared  operator 
in  the  numerator  yields 

m  n _ 

m  n 

where  „ 

In  order  to  simplify  the  expression,  we  will  now  use  the  approximation  that  the 
values  ^  are  a  constant  for  all  combinations  of  m  and  n.  This  approximation  is 
somewhat  justified  for  the  following  two  reasons.  One  is  that  the  phase  of  the  ’s  has 
been  canceled  due  to  the  magnitude  operators.  The  other  is  that  we  have  tended  to 
examine  average  matching  coefficients  for  a  grid  of  source  locations  in  a  given  region. 
This  tends  to  average  out  the  effects  seen  because  some  modes  are  excited  to  different 
extents  due  to  the  placement  of  the  source  in  depth.  Further  ramifications  of  this 
assumption  will  be  discussed  shortly.  Simplifying  Eq.  3.9  with  the  assumption  just  stated 
gives  us 


_  m _ rj_ 


(3.9) 


m  n 


(3.10) 


where  M  is  the  number  of  modes  and  we  will  call  Yg  an  “obliquity  characterization 
number.”  Ideally  this  measure  would  give  us  the  expected  matching  coefficient  for  a  given 
environment.  Note  that  if  all  the  0^’s  are  equal  to  one,  effectively  giving  us  a  perfect 
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Table  3.3.  Matching  Coefficients  V5.  Obliquity  Characterization  Numbers. 


Example  # 

Obliquity 
Characterization 
Number  (Yg) 

Matching 

Coefficient 

3C 

0.999 

0.986 

3A 

0.997 

0.984 

3D 

0.994 

0.979 

3E 

0.983 

0.966 

3B 

0.965 

0.944 

extrapolation,  that  Tg  is  also  equal  to  one  which  is  the  maximum  value  of  the  matching 
coefficient. 

We  will  now  experimentally  explore  the  relationship  between  the  obliquity 
characterization  number  and  the  matching  coefficient.  Note  that  we  are  still  assuming 
perfect  knowledge  of  the  environmental  parameters,  as  reflected  in  our  modeled  fields.  In 
addition  to  the  two  examples  already  run,  three  more  environmental  scenarios  were 
devised  in  order  to  vary  the  obliquity  characterization  number  (as  well  as  the  matching 
coefficient).  Examples  3C,  3D,  and  3E  utilized  the  same  general  geometry  as  Exs.  3A  and 
3B,  but  the  constant  sound  speed  was  changed  to  1465  m/s,  1535  m/s,  and  1560  m/s, 
respectively.  Table  3.3  contains  the  obliquity  characterization  numbers  and  average 
matching  coefficients  for  all  the  examples.  Note  that  the  table  entries  are  listed  in  order  of 
decreasing  Yg .  Notice  the  consistent  decreasing  nature  of  the  matching  coefficient  as  the 
obliquity  characterization  number  decreases. 

Figure  3.9  shows  the  data  points  from  Table  3.3  along  with  a  least-squares  linear 
fit.  The  equation  for  this  least-squares  fit  is 

1.236 Yg-0.248. 


(3.11) 
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Figure  3.9.  Matching  Coejf.  V5.  Obliquity  Char.  Number  and  Linear  Fit. 

As  can  be  seen  from  examination  of  Fig.  3.9,  the  linear  fit  is  an  excellent  one  to  the  data 
points.  Note  that  extrapolating  the  quadratic  fit  out  to  an  obliquity  characterization  number 
of  one  results  in  a  matching  coefficient  of  approximately  only  0.988.  Also,  the  ideal  slope 
of  the  linear  fit  should  be  one. 

The  fact  that  the  extrapolated  line  does  not  produce  a  matching  coefficient  of  one 
for  an  obliquity  characterization  number  of  one  can  be  explained  by  the  numerical 
integration  method  used  to  do  the  extrapolation.  This  will  be  explored  in  Ch.  5,  and  the 
relationship  between  Yg  and  l/r|  will  be  revisited  there.  The  variation  of  the  slope  of  the 
linear  fit  from  one  is  due  to  the  assumption  that  the  values  of  „  were  constant  for  all  m 
and  n.  These  values  are,  in  fact,  not  constant  and  tend  to  decrease  as  m  and  n  increase  due 
to  the  factor  in  the  representation  of  the  pressure  field.  This  consistent 

nonuniformity  of  the  ^’s  is  then  the  cause  of  the  steeper  slope  in  the  fitted  line.  This 
effect  could  be  incorporated  into  a  new  obliquity  characterization  number  calculation,  but 
this  would  involve  a  double  summation  in  the  denominator  and  would  be  less  intuitive  to 
analyze.  As  it  is,  the  simply  stated  expression  for  the  obliquity  characterization  number 
gives  a  good  relative  measure  of  how  the  matching  coefficient  in  an  extrapolation  is 
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affected  by  the  obliquity  factors. 

3.5.  Obliquity  Factor  Effects:  Attenuating  Sediment  Layer  (Example  3F) 

A  common  situation  for  shallow  water  acoustic  propagation  is  to  have  an 
attenuating  sediment  layer  beneath  the  water  column,  as  opposed  to  the  rigid  boundary 
that  has  been  assumed  for  the  previous  examples.  We  will  use  the  environmental 
parameters  depicted  in  Fig.  3.10.  The  sediment  parameters  were  determined  from 
experimental  data  taken  from  a  deep  ocean  basin  off  the  west  coast  of  Vancouver  Island, 
B.C.,  and  are  representative  of  a  sediment  with  a  “clayey-silt”  composition  [Dosso  and 
Chapman].  The  somewhat  arbitrary  negative  gradient  in  the  water  column’s  sound  speed 
profile  produces  a  downward  refracting  effect  on  the  sound  traveling  through  the  water,  in 
effect  “forcing”  the  sound  into  more  interaction  with  the  sediment. 

The  array  geometry  for  Ex.  3F  will  be  similar  to  that  for  Exs.  3A  and  3B  in  that  the 
calibration  range  will  be  =  20  km ,  and  the  arrays  will  be  populated  with  ideal  point 
sources  with  spacing  AZj  =  =  1  m . 

3.5.1.  Modal  Structure 

Using  a  source  frequency  of  300  Hz,  KRAKEN  calculates  52  propagating  modes 
for  the  scenario  depicted  in  Fig.  3.10.  Table  3.4  contains  a  subset  of  the  horizontal 
wavenumbers  (eigenvalues)  and  associated  approximate  grazing  angles  for  this  scenario. 
Since  attenuation  due  to  material  absorption  can  be  dealt  with  by  allowing  the  sound  speed 
[Jensen  et  al],  or,  equivalently,  the  wavenumber  [Boyles],  to  become  complex,  the 
eigenvalues  for  this  scenario  are  complex. 

The  mode  functions  corresponding  to  the  eigenvalues  listed  in  Table  3.4  are  shown 
in  Fig.  3.11.  Note  that  the  sediment  layer  begins  at  a  depth  of  150  m,  and  that  significant 
extension  of  the  mode  functions  into  the  sediment  begins  somewhere  between  modes  15 


and  22. 
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Figure  3.10.  Environment  with  attenuating  sedimentary  layer. 
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Table  3.4.  Eigenvalues  for  Ex.  3F. 


Mode  #  (m) 

Eigenvalue  (^^) 

Grazing  Angle  (xj 
(deg.) 

2 

1.267  +  19.201x10-^ 

5.8 

4 

1.262  +  18.846x10-^ 

7.6 

6 

1.259  +  19.218x10-^ 

8.8 

8 

1.254 +  i  1.338x10-^ 

10.0 

10 

1.248  +  1  1.967x10-^ 

11.4 

15 

1.228  +  15.303x10-^ 

15.5 

22 

1.193 +  i  1.306x10-^ 

20.5 

29 

1.162  +  11.440x10-^ 

24.2 

36 

1.131  +i  1.563x10-^ 

27.4 

43 

1.101  +i  1.466x10-^ 

30.2 

Referring  back  to  the  normal  mode  equation  governing  range-independent 
propagation  (Eq.  1.10),  we  see  that  the  contribution  of  each  mode  to  the  acoustic  pressure 
at  a  given  range  r  is  weighted  by  an  e  term.  When  the  eigenvalues  are  complex,  i.e., 
+  ib^  (b^  is  defined  to  be  a  positive  number),  the  weighting  factor  becomes 


+  ia^r  -b^r  ia^r 

e  =  e  =  e  e  -  (x{r)e 


(3.12) 


where  a(r)  is  a  multiplicative  factor  that  decreases  with  increasing  range,  i.e.,  the 
attenuation  due  to  this  factor  increases  with  range.  Also,  this  factor  decreases  as  b^^ 
increases.  In  Table  3.4,  we  see  that  the  imaginary  part  of  the  eigenvalues  increases  by 
roughly  two  orders  of  magnitude  somewhere  between  mode  numbers  15  and  22.  This 
makes  intuitive  sense  in  that  this  is  also  the  region  where  the  mode  functions  begin  to 
extend  into  the  sediment,  i.e.,  modes  that  interact  more  with  the  sediment  are  more 


Depth  (m)  Depth 
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Figure  3.12.  Attenuation  at  20  km  due  to  imaginary  parts  of  eigenvalues  (sediment). 

severely  attenuated  in  range.  Figure  3.12  shows  the  attenuation  due  solely  to  the  imaginary 
parts  of  the  eigenvalues  in  this  scenario  at  a  range  of  20  km.  We  see  attenuations  of  greater 
than  60  dB  for  all  modes  after  number  18,  and  attenuations  greater  than  120  dB  for  all 
modes  after  number  24. 

Keeping  these  issues  in  mind,  Fig.  3.13  shows  a  transmission  loss  plot  for  our 
sediment  interaction  scenario.  Note  the  high  loss  in  the  sediment  (greater  than  90  dB)  for 
ranges  greater  than  5  km.  The  modes  that  have  the  most  penetration  into  the  sediment  have 
been  “stripped”  by  this  point  in  range,  and  the  modes  that  have  little  or  no  penetration  into 
the  sediment  continue  to  propagate  (there  is  some  incidental  acoustic  energy  propagating 
in  the  sediment  for  ranges  >  5  km). 

3.5.2.  Averaged  Extrapolation 

We  can  now  examine  what  effect  mode  stripping  by  the  sediment  layer  has  on  the 
acoustic  field  extrapolation  algorithm.  As  was  seen  in  Ex.  3B,  the  existence  of  very  steeply 
propagating  modes  (leading  to  extreme  obliquity  factors)  caused  a  degradation  of  the 
extrapolation  algorithm.  One  would  expect,  then,  that  the  elimination  of  steeply 
propagating  modes  would  aid  the  extrapolation  algorithm,  and  this  is  indeed  the  case. 


Figure  3. 
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In  order  to  find  a  representative  matching  coefficient  for  this  scenario,  source 
positions  were  varied  over  a  grid  of  ranges  from  1  km  to  5  km  spaced  at  1  km  and  depths 
from  15  m  to  135  m  spaced  at  15  m.  The  average  value  of  the  matching  coefficient  over 
this  grid  of  source  locations  was  approximately  0.996  with  a  variance  of  1.764x10'^. 
Figure  3.14  shows  the  averaged  extrapolation  results  for  a  hypothetical  source  location 
=  (5  km,  30  m)  .  The  matching  coefficient  is  approximately  0.996.  The 
primary  and  alternate  extrapolations  are  not  shown  because  the  obliquity  factors  have  a 
negligible  effect,  making  the  results  of  these  two  extrapolations  nearly  identical. 

3.5.3.  Obliquity  Characterization  of  the  Medium 

Example  3F  is  significantly  different  from  Examples  3A  and  3B  in  that  not  all  of 
the  calculated  propagating  modes  contribute  significantly  to  the  fields  at  the  receiving 
array.  Based  on  Fig.  3.12,  it  seems  reasonable  to  assume  that  the  first  18  modes  contribute 
the  bulk  of  the  acoustic  energy  to  the  field  measured  along  the  receiving  array,  so  we  will 
use  only  these  modes  in  the  calculation  of  the  obliquity  characterization  number,  which 
turns  out  to  be  approximately  1.0.  This  is  the  expected  result  due  to  the  shallow  grazing 
angles  of  the  propagating  modes,  resulting  in  obliquity  factors  very  close  to  one.  This  in 
turn  helps  to  explain  the  high  matching  coefficient  between  the  ground  truth  and 
extrapolated  fields  of  Ex.  3F. 

Also  note  that  there  is  not  as  great  a  discrepancy  between  the  predicted  matching 
coefficient  for  an  obliquity  characterization  number  of  one.  In  Sec.  3.4,  the  linear  fit  was 
extrapolated  to  a  value  of  j.K’l  =  0.988  for  Yg  =  1 ,  whereas  in  this  example  an  actual 
value  of  0.996  was  obtained  for  Yg  =  1 .  This  can  again  be  explained  by  the  numerical 
integration  method  used,  and  will  be  revisited  in  Ch.  5. 

3.6.  Summary  of  Chapter  3 

Chapter  3  was  intended  mainly  as  an  illustration  through  numerical  modeling  of 
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the  mathematical  relationships  explored  in  Ch.  2.  Through  qualitative  comparisons  it  was 
seen  how  the  overemphasis  of  the  individual  modes  in  the  primary  extrapolation  and 
underemphasis  in  the  alternate  extrapolation  lead  to  distinctly  different  types  of  distortions 
in  the  extrapolated  fields.  It  was  also  seen  how  averaging  these  two  estimates  ameliorates 
these  effects  and  produces  a  better  approximation  of  the  ground  truth  field.  It  was  shown 
how  a  quantitative  comparison  of  the  extrapolated  and  ground  truth  fields  can  be  made 
through  use  of  the  matching  coefficient,  and  how  approximations  of  the  obliquity  factors 
for  each  mode  in  a  given  environmental  scenario  can  give  an  approximation  to  the 
expected  matching  coefficient. 

Chapter  4  will  examine  the  effects  of  inaccurate  knowledge  of  the  environmental 
parameters  used  to  produce  the  modeled  fields  on  the  extrapolation  algorithm. 


Chapter  4:  Acoustic  Field  Extrapolation  -  Modeled  Field  Mismatch 

4.1.  Introduction  to  Chapter  4 

The  examples  given  in  Ch.  3  assumed  perfect  knowledge  of  the  medium  to  be  used 
as  input  to  an  acoustic  field  modeling  algorithm  in  order  to  calculate  p .  This  led,  then,  to 
“exact”  modeled  fields  (in  the  sense  that  they  reflect  what  the  actual  fields  would  look  like 
along  the  reference  array),  and  extrapolated  fields  that  are  affected  only  by  the 
extrapolation  algorithm  itself.  These  effects  were  seen  to  be  modal  dependent  obliquity 
factors,  which  became  more  severe  for  more  steeply  propagating  energy. 

We  will  now  examine  the  effects  of  using  erroneous  environmental  information  in 
the  modeling  step,  leading  to  a  situation  we  will  refer  to  as  “environmental  mismatch.” 
The  erroneous  environmental  information  results  in  modeled  fields  that  do  not  accurately 
reflect  what  the  actual  field  would  look  like  along  the  reference  array,  and  therefore 
degrade  the  extrapolation  results.  Numerous  studies  of  the  effect  of  mismatch  on  matched- 
field  processing  algorithms  have  been  done  (see  references  in  [Baggeroer  et  al],  as  well  as 
[Gingras]).  This  chapter  will  examine  the  effect  of  mismatch  due  to  errors  in  the  sound 
speed  profile  of  the  water  column  as  well  as  errors  in  sediment  parameters. 

This  chapter  will  take  a  “perturbational”  approach  to  the  problem,  i.e.,  for  a  small 
perturbation  in  the  water  column  sound  speed  profile,  what  is  the  change  in  the  expression 
for  the  modeled — and  subsequently  the  extrapolated — field?  This  will  be  accomplished 
with  the  normal  mode  theory  used  throughout  this  work.  First  an  existing  derivation  from 
the  field  of  quantum  mechanics  will  be  presented  that  shows  how  an  equation  like  the 
modal  equation  can  be  treated  perturbationally.  This  result  will  then  be  translated  into  the 
underwater  acoustic  modal  equation  form  in  order  to  show  the  changes  in  the  horizontal 
wavenumbers  and  mode  functions  due  to  perturbations  in  the  sound  speed  profile  in  the 
water  column.  A  “zero-order”  acoustic  waveguide  scenario  will  be  used  to  determine 
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some  general  characteristics  of  the  first-order  perturbations.  The  perturbations  will  be 
integrated  into  the  field  extrapolation  expressions,  and  examples  of  errors  in  the  water 
column  sound  speed  profile  as  well  as  the  sediment  properties  wilt  be  given  and  their 
effects  on  the  algorithm  shown. 

4.2.  Quantum  Mechanical  Perturbation  Theory 

In  quantum  mechanics,  the  time-independent  Schrodinger  equation  has  a  form 
very  similar  to  our  modal  equation  for  underwater  sound  propagation.  The  consequences 
of  perturbing  the  “energy  operator”  in  quantum  mechanics  are  well  known,  and  can  serve 
as  a  guide  in  doing  a  similar  operation  with  the  modal  equation.  The  following  first-order 
perturbation  derivation  follows  the  order  and  notation  found  in  [Yariv].  A  parallel 
development  can  be  found  in  [Schiff]. 

The  time-independent  Schrodinger  equation  is  given  as 


Hnu  =  E  u  , 

^  m  m 


(4.1) 


where  the  are  the  wavefunctions  (eigenfunctions)  for  a  given  state,  and  the  are  the 
energies  (eigenvalues)  for  that  state.  Hq  is  the  energy  operator,  or  Hamiltonian  (the  exact 
form  of  the  Hamiltonian  is  not  of  immediate  relevance  here).  The  problem  posed  is  thus:  if 
the  Hamiltonian  is  perturbed  from  Hq  to  Hq  +  H' ,  what  are  the  new  eigenfunctions  and 
eigenvalues?  Assuming  Hq  »  H' ,  we  can  use  the  following  approach. 

We  set  the  perturbed  Hamiltonian  to  be  Hq  +  XH' ,  where  X  is  simply  a  parameter 
that  is  set  to  one  where  a  perturbation  is  desired  and  set  to  zero  where  it  is  not.  For  the  final 
results,  X  will  be  set  to  one.  We  are  now  looking  for  the  eigenvalues  W  and  eigenfunctions 
4^  that  satisfy 


(Hq  +  XH')'¥  =  WT. 


(4.2) 
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Expanding  W  and  'F  in  power  series  in  X,  i.e., 


'F  =  'Fo  +  ^.'Fj  +  X  'F2+"- 


W  =  Wq  +  XWj  +X  + 


and  substituting  these  expansions  back  into  Eq.  4.2,  we  obtain 


{Hq  +  XH')  ('Fo  +  ?i'Fj+?i^'F2+...)  =  (Wo  +  >.W,  +  ?.V2  +  ...)  X 

('Fo  +  ?.'Fj+?i^'F2 +...). 


Now  we  can  equate  like  powers  of  X  on  both  sides  of  Eq.  4.4,  giving 


and  so  on.  Comparing  the  first  line  of  Eq.  4.5  with  Eq.  4.1,  we  see  that  the  zero-order 
solutions  are 

^0  =  ^m 

^0  =  (4.6) 

or  the  unperturbed  eigenfunctions  and  eigenvalues. 

To  determine  the  first  order  perturbations  to  the  eigenfunctions  and  eigenvalues, 
'Fj  and  Wj  respectively,  we  first  expand  'Fj  in  terms  of  the  as 


'i'l  = 


(4.7) 


where  the  are  coefficients  for  the  first-order  perturbation.  Noting  that  can  be 
substituted  for  Hq  (see  Eq.  4.1)  and  using  Eq.  4.7,  the  second  line  of  Eq.  4.5  (A,  ,  first- 
order  perturbation)  can  now  be  written 

+H'u  =E  +W,u  .  (4.8) 

^  n  n  n  m  mZmj  n  n  \  m  ^  ^ 

n  n 

The  eigenfunctions  here  take  on  an  orthonormality  property  similar  to  the  mode 
functions  for  the  underwater  acoustic  propagation  formulas,  i.e., 

[m  M  =  6  ,  (4.9) 

]  m  n  mn^ 

where  the  *  denotes  conjugation  and  6^^  is  a  Kronecker  delta  function.  Taking  advantage 

* 

of  this,  we  premultiply  Eq.  4.8  by  and  integrate,  resulting  in 

+  (»;,  H’uJ  =  (4.10) 

where  {a,  b)  denotes  an  inner  product  of  a  and  b.  Letting  k  =  mm  Eq.  4.10  yields 


H'u  ) 


(4.11) 


giving  us  the  first  order  correction  for  the  energy  E  .  For  k^m,'we  find  that 


(m^,  H'uJ 


k^m . 


(4.12) 


It  can  be  shown  [Yariv]  that  =0,  so  the  complete  first-order  perturbational 
expressions  for  the  eigenfunctions  and  eigenvalues  are  given  by 
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4.3.  Underwater  Acoustic  Modal  Equation  Perturbation  Theory 

We  can  now  apply  the  above  results  to  the  underwater  acoustic  modal  equation  to 
determine  what  the  changes  in  the  mode  functions  and  horizontal  wavenumbers  would  be 
for  perturbations  in  the  water  column  sound  speed  profile.  Similar  results  are  shown  in 
[Rajan  et  al],  and  we  use  elements  of  their  notation  here.  A  parallel  derivation  was  also 
done  in  [Rouseff  (Feb.  1993)]. 

We  begin  by  assuming  a  piecewise  constant  density  profile,  and  rewriting  the 
modal  equation,  Eq.  1.4,  as 


where 


(4.15) 


is  the  unperturbed  depth-dependent  acoustic  wavenumber,  and  Cq(z)  is  the  unperturbed 
sound  speed  profile  in  the  water  column.  Comparing  Eq.  4.14  with  Eq.  4.1,  we  see  that  the 
mode  functions  correspond  to  the  quantum  mechanical  wavefunctions 

(eigenfunctions)  u^,  the  squared  horizontal  wavenumbers  (^^  )  correspond  to  the 
quantum  mechanical  energy  levels  (eigenvalues)  E^,  and  the  operator  term  in  [  ]’s 
corresponds  to  the  Hamiltonian  Hq.  We  will  use  the  same  notation  and  terminology  for 
this  operator  as  was  used  in  the  quantum  mechanical  derivation,  although  the  two 
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operators  have  different  forms  when  expanded. 

Since  it  appears  explicitly  in  the  perturbed  eigenfunction  and  eigenvalue 
expression,  the  next  task  is  to  determine  the  perturbation  in  the  Hamiltonian,  H' ,  due  to  a 
perturbation  in  the  sound  speed  profile.  We  will  denote  the  sound  speed  perturbation  as 
Ac(z) ,  and  assume  that  Ac(z) «  Cq(z)  .  We  now  write  the  expression  for  the  perturbed 
depth-dependent  acoustic  wavenumber  as 


k{z)  = 


0)  _  _  CO  _  coAc(z) 

Cq{z)  +  Ac(z)  ^2^^^  _  (Ac(z))  ^  Az) 


(4.16) 


assuming  terms  of  (Ac(z))  are  negligible.  We  then  find  that 


;  2/  X  _  0) 

k  (z)  =  -j— 


2(0  Ac(z) 


(O^Ac(z))" 


(z))  -2{k' 


^Ac(z) 


(4.17) 


2 

again  neglecting  terms  of  ( Ac(z))  .  Inserting  this  into  the  expression  for  the  Hamiltonian, 
we  see  that 


Ho  +  H'  =  -^  +  k\z)  =  -^+ 
dz^  dz" 


(4.18) 


Comparing  this  to  Eq.  4.14,  we  can  see  that  the  perturbation  in  the  Hamiltonian  is 


H'  =  -2(k^^\z))^^, 

Cq{z) 


(4.19) 


Note  also  from  Eq.  4.14  that  the  energies  (eigenvalues)  of  the  quantum  mechanical 
expression  are  the  counterpart  of  the  square  of  the  horizontal  wavenumbers  for  the 
underwater  acoustic  modal  equation.  We  can  say  that 


75 


=  )  =(^m  +  )  "  (S;;,  ) 


(4.20) 


neglecting  terms  of  ( A^^*  ^ )  ,  and  consequently 


A4r>  - 


24 


(0) 


(4.21) 


Translating  the  quantum  mechanical  perturbation  results  (Eq.  4.10),  we  see  to  first 
order  that 


A4i  =  (vrw.*>rft))'  (4.22) 

Keeping  in  mind  that  the  inner  product  orthonormality  condition  for  the  modal  equation 
involves  a  weighting  by  the  inverse  of  the  density  profile  (see  Eq.  1.6),  we  can  combine 
Eqs.  4.21  and  4.22  to  determine  the  first  order  perturbation  expression  for  the  horizontal 
wavenumbers  and  write  it  as 


0  ° 


(4.23) 


Translating  the  proper  quantum  mechanical  results,  we  can  also  write  the 
expression  for  the  perturbation  in  the  mode  functions  as 


(D.JOK 


(4.24) 
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where 


a 


(1) 

n 


(yf ’fe). 


n^m 


0 


n  =  m. 


Substituting  in  the  proper  expressions  results  in 


(4.25) 


-2W„(Z) 


=  I 


n^m 


J  P”’  ao  (k  iz')) 

0 _ 

iC)'-(C)' 


(4.26) 


4.4.  Discussion  of  Perturbation  Results 

We  will  now  examine  some  of  the  implications  of  the  above  perturbation  results. 
We  will  use  a  “zero-order”  acoustic  waveguide  case  in  order  to  gain  some  insight  into  how 
perturbations  in  the  sound  speed  profile  affect  the  horizontal  wavenumbers  and  mode 
functions.  Both  constant  offsets  and  random  fluctuations  will  be  examined. 


4.4.1.  Constant  Offset 

Examples  3A  and  4A  will  serve  as  our  “zero-order”  acoustic  waveguide  examples. 
In  Ex.  3A  we  had  an  idealized  ocean  with  a  pressure  release  surface  and  a  perfectly  rigid 
bottom  at  a  depth  of  h  =  50  m ,  and  we  used  a  constant  sound  speed  profile  in  the  water 
column  of  1500  m/s.  Note  that  the  mode  functions  and  horizontal  wavenumbers  can  be 
calculated  analytically  for  this  case  (see  Eq.  3.1).  Example  4A  will  use  a  similar  scenario 
with  a  constant  sound  speed  of  1510  m/s. 

We  can  examine  our  first-order  perturbational  results  in  order  to  see  what  the 
effects  of  perturbing  the  sound  speed  profile  (as  in  going  from  Ex.  3A  to  Ex.  4A)  on  the 
horizontal  wavenumbers  and  mode  functions  should  be.  Examining  Eq.  4.23,  we  realize 
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that  the  density,  depth-dependent  wavenumber,  and  sound  speeds  (both  the  Ac(z)  and  the 
Cq(z))  are  constants  with  depth  which  can  be  pulled  out  of  the  integral.  Taking  the 
orthonormality  condition  into  account,  we  see  that  the  perturbation  in  each  horizontal 
wavenumber  should  be 


where  A  is  a  constant  for  all  m  (A  is  positive  for  Ac(z)  >  0).  In  other  words,  the  change  in 
each  horizontal  wavenumber  is  a  constant  divided  by  the  unperturbed  horizontal 
wavenumber.  Since  the  eigenvalues  are  ordered  such  that  , ,  the  perturbations 

should  be  greater  as  m  increases.  This  is  borne  out  in  Table  4.1  which  lists  the  horizontal 
wavenumbers  (to  four  decimal  places)  for  Examples  3A  and  4A  along  with  their 
differences.  The  perturbations  do  indeed  grow  as  m  increases. 

Using  the  first  mode  in  order  to  get  an  estimate  of  the  constant  A  (without  rounding 
the  KRAKEN  output),  we  see  that 

A  =  =2.61449  X  10“^  (4.28) 

Using  this  to  predict  what  should  be,  we  arrive  at  an  approximate  value  of -0.0133. 

Keeping  in  mind  that  this  is  a  first  order  approximation,  it  is  quite  close  to  the  actual  value 
of -0.0138. 

We  can  do  a  similar  simplification  for  the  mode  function  perturbations  of  Eq.  4.26. 
This  time,  however,  when  the  constant  terms  are  pulled  outside  the  integral,  we  are  left 
with  inner  products  of  orthogonal  components  (since  the  sum  is  over  all  n^m ),  and  the 
result  is  that  A\|/^'\z)  =  0.  This  coincides  with  the  result  from  Ch.  3  where  it  was  seen 
that  the  mode  functions  are  indeed  the  same  for  simple  constant  offsets  of  the  sound  speed 
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Table  4.1.  Eigenvalue  Comparison  for  Ex.  4A. 


Mode  #  (m) 

Example  3A 

Example  4A 

Difference 

1 

.6275 

.6234 

-.0041 

2 

.6212 

.6170 

-.0042 

3 

.6084 

.6041 

-.0043 

4 

.5886 

.5841 

-.0045 

5 

.5611 

.5564 

-.0047 

1 

6 

.5247 

.5198 

-.0049 

7 

.4775 

.4720 

-.0055 

8 

.4156 

.4093 

-.0063 

9 

.3310 

.3230 

-.0080 

10 

.1962 

.1824 

-.0138 

for  environments  of  this  type. 

4.4.2.  Random  Fluctuations 

We  will  now  examine  the  situation  where  the  perturbation  to  the  sound  speed  is  a 

random  process.  Specifically,  Ac(z)  will  be  a  wide-sense  stationary  zero-mean  white 

2 

Gaussian  process  with  variance  a  .  For  this  case,  we  can  write  the  first-order 
perturbations  to  the  horizontal  wavenumbers  as 

0 


Taking  the  expected  value  [Papoulis]  of  this  expression  yields 
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1  '  ^ 

=  =§;ji'vT(z'))  E{Ac(z'))dz' 

=  0  (4.30) 

since  Ac(z)  is  zero-mean. 

We  can  also  calculate  the  variance  of  the  perturbations  of  the  horizontal 
wavenumbers  as 
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Table  4.2.  Results  of  Random  Fluctuations  on  Eigenvalues. 


Mode  #  (m) 

Unperturbed 

Sample  Mean 

Sample  Variance 

1 

.628 

.628 

5.23x10'^ 

2 

.621 

.621 

5.35x10'^ 

3 

.608 

.608 

5.76x10'^ 

4 

.589 

.589 

6.53x10'^ 

5 

.561 

.561 

6.41x10'^ 

6 

.525 

.525 

7.46x10'^ 

7 

.477 

.477 

8.86x10'^ 

8 

.416 

.416 

1.29x10'^ 

9 

.331 

.331 

1.83x10'^ 

10 

.196 

.196 

5.45x10*^ 

since  it  can  be  shown  that  the  integral  of  the  mode  functions  to  the  fourth  power  is  a 

constant  for  all  m,  given  the  form  of  the  mode  functions  in  Eq.  3.1. 

While  the  mean  of  the  perturbations  of  the  horizontal  wavenumbers  is  zero  for  all 

modes,  the  variance  grows  proportionally  to  the  inverse  of  the  square  of  the  unperturbed 

horizontal  wavenumbers.  A  Monte  Carlo  experiment  was  run  using  KRAKEN  to  calculate 

the  horizontal  wavenumbers  for  realizations  of  perturbed  sound  speed  profiles. 

2 

Realizations  of  Ac(z)  (spacing  in  z  was  1  m)  with  a  =  1  m/s  were  added  to  a  constant 
sound  speed  profile  of  1500  m/s  (i.e.,  perturbing  the  environmental  parameters  of  Example 
1).  The  results  over  500  realizations  are  listed  in  Table  4.2.  Note  that  the  sample  means 
match  the  unperturbed  horizontal  wavenumbers  for  all  modes,  and  the  variances  show  a 
more-or-less  consistent  increase  with  m  (there  is  an  anomaly  at  m  =  5). 

Doing  a  similar  operation  as  was  done  for  determining  the  constant  A,  we  can  use 
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the  results  from  the  first  mode  to  determine  A" \ 

A"- (^[°^)^s2.06xl0"^  (4.34) 

Using  this  to  predict  what  the  variance  should  be  for  the  tenth  mode  results  in  a  value  of 
5.37x10'^,  which  compares  very  favorably  with  the  actual  value  of  5.45x10'^.  Again  we 
see  the  trend  in  which  the  higher  order  modes  are  affected  to  a  greater  degree  by 
perturbations  in  the  sound  speed  profile.  The  expression  for  the  variance  of  the  A\|/^  ^  (z) 
is  easily  derivable,  but  is  of  little  use  for  gaining  insight  into  how  the  modes  are  affected 
by  the  sound  speed  fluctuations. 

It  is  important  to  note  that  the  results  in  this  section  (4.4)  apply  only  to  the  special 
“zero-order”  acoustic  waveguide  case,  whereas  the  general  results  from  Sec.  4.3  apply  to 
any  environmental  scenario  that  satisfies  the  assumptions  stated  there.  By  making  the 
background  sound  speed  constant  with  depth,  we  were  able  to  move  terms  based  on  it 
outside  the  depth  integral,  and  some  illustrative  analytical  results  were  obtained.  For  other 
scenarios,  Eqs.  4.23  and  4.26  could  be  numerically  integrated  in  order  to  provide  estimates 
of  the  changes  in  the  horizontal  wavenumbers  and  mode  functions  for  a  given  perturbation 
of  the  sound  speed. 

4.5.  Mismatch  Effects:  Water  Column  Sound  Speed  Profile  (Example  4B) 

This  section  will  provide  an  example  of  the  effects  on  the  acoustic  field 
extrapolation  algorithm  due  to  mismatch  in  the  sound  speed  profile  of  the  water  column. 
Random  fluctuations  about  the  true  sound  speed  profile  will  be  considered  here. 

4.5.1.  Environmental  Parameters  and  Modal  Structure 

For  this  example,  we  will  assume  imperfect  knowledge  of  the  environment  in  the 
form  of  a  mismatched  sound  speed  profile  in  the  water  column.  The  “true”  environmental 


82 


Figure  4.1.  "True”  and  perturbed  sound  speed  profile  in  water  column. 

conditions  will  be  the  same  as  those  found  in  Sec.  3.4,  i.e.,  an  attenuating  sediment  layer 

underlying  the  water  column.  We  choose  this  environment  in  order  to  focus  attention  on 

the  effects  due  to  mismatch  instead  of  the  obliquity  factors.  The  “calibration”  data,  then, 

will  be  taken  (simulated,  in  this  case)  using  these  conditions.  The  erroneous  environmental 

conditions  we  will  use  in  constructing  our  modeled  fields  will  be  the  same  as  the  “true” 

conditions  except  for  random  fluctuations  Ac(z)  in  the  sound  speed  of  the  water  column. 

Ac(z)  will  be  a  wide-sense  stationary  zero-mean  white  Gaussian  process  with  variance 
2 

a  =  1  m/s .  A  plot  of  the  true  and  perturbed  sound  speed  profiles  can  be  seen  in  Fig.  4.1. 

This  mis-estimation  of  the  sound  speed  profile  leads  to  perturbations  in  the 
modeled  mode  functions  and  horizontal  wavenumbers  as  was  seen  in  the  previous  section. 
Although  the  mean  values  of  the  perturbations  were  shown  to  be  zero,  any  given 
realization  of  Ac(z)  will  impart  a  finite  perturbation  to  the  mode  functions  and  horizontal 
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Table  4.3.  Eigenvalue  Comparison  for  Ex.  4B. 


Mode  #  (m) 

Eigenvalue 

Change  in  Eigenvalue 

2 

1.267  +  i  9.475x10-^ 

-1.945x10-'^  +  i  2.739x10-'^ 

4 

1.262  +  17.404x10-^ 

1.645x10-^-1 1.442x10-^ 

6 

1.259  +  19.951x10-^ 

-6.378x10-^  +  17.328x10-^ 

8 

1.254 +  i  1.489x10-^ 

3.672x10-^  +  1  1.515x10-^ 

10 

1.248 +  i  1.988x10-^ 

-9.394x10-^  +  i  2.017x10-'^ 

15 

1.228  +  15.270x10-^ 

5.519x10-3-13.296x10-'^ 

22 

1.193  +i  1.299x10-^ 

2.741x10-^-16.843x10-^ 

29 

1.162 +  i  1.421x10-3 

-2.694x10-3  _  •  1.908x10-3 

36 

1.131  +i  1.570x10-3 

-2.742x10-^  +  i  6.573x10-^ 

43 

1.101  +i  1.470x10-3 

-5.245x10-^  +  i  3.560x10-^ 

wavenumbers.  Table  4.3  shows  a  comparison  between  the  “true”  and  perturbed  horizontal 
wavenumbers  for  this  example.  Although  it  is  difficult  to  draw  any  general  conclusions,  it 
does  appear  that  greater  variations  in  the  horizontal  wavenumbers  appear  in  the  lower 
order  modes. 

This  can  be  explained  in  the  following  manner.  Referring  back  to  Eq.  4.23,  we  see 
that  the  perturbation  for  a  given  horizontal  wavenumber  is  expressed  as  an  integral  over 
the  squared  mode  function  times  environmental  parameters,  including  the  sound  speed 
perturbations.  Since  the  Ac(z)  are  zero  for  depths  past  the  water/sediment  boundary,  the 
integral  need  only  go  from  0  to  h  =  150m.  Also  keep  in  mind  that  due  to  the 
orthonormality  condition,  each  mode  has  unit  “energy”  when  their  squares  are  integrated 
from  0  to  «>  (or  in  this  case,  to  the  attenuating  basement  at  h  =  300  m ).  Since  higher 
order  mode  functions  have  increasingly  more  “energy”  in  the  sediment,  and,  hence,  less 
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“energy”  that  is  encompassed  by  the  integral  of  Eq.  4.23,  we  would  expect  that  the  lower 
order  modes  would  be  affected  more  by  the  water  column  sound  speed  fluctuations  due  to 
the  integral  part  of  Eq.  4.23.  This  effect  ends  up  competing  with  the  inverse  weighting  of 
the  unperturbed  horizontal  wavenumbers,  but  since  for  this  case  the  horizontal 
wavenumbers  vary  over  a  relatively  small  range,  the  integration  factor  wins  out. 

Figure  4.2  shows  a  plot  comparing  the  mode  functions  for  the  “true”  and  perturbed 
cases  (the  perturbed  modes  are  plotted  in  gray).  Small  errors  can  be  observed  in  the  first 
four  mode  functions  shown,  but  the  rest  show  near  exact  coincidence  between  the  “true” 
and  perturbed  mode  functions.  Again  we  can  appeal  to  our  perturbation  relationships. 
Equation  4.26  shows  the  first-order  perturbation  results  for  the  mode  functions.  Note  that 
for  a  given  mode  m,  the  perturbation  to  the  mode  function  is  a  sum  of  integrals  (again  only 
over  the  water  depth)  over  that  weighted  mode  function  times  the  rest  of  the  mode 
functions  to  which  it  is  orthogonal  (the  weighting  skews  the  orthogonality).  Since  the 
higher  order  modes  have  less  “energy”  in  the  water  column,  we  can  expect  that  they  will 
be  perturbed  less  than  the  lower  order  modes  which  have  nearly  all  their  “energy”  in  the 
water  column. 

4.5.2.  Extrapolation  Results 

An  extrapolation  experiment  was  run  using  the  same  array  geometry  as  in  Ex.  3F. 
Extrapolations  were  performed  for  a  grid  of  source  locations  varying  in  range  from  0.2  km 
to  10  km  spaced  at  0.2  km  and  in  depth  from  15  m  to  135  m  spaced  at  5  m.  Instead  of 
averaging  over  the  entire  grid,  the  matching  coefficients  for  a  given  range  were  averaged, 
yielding  an  average  matching  coefficient  vs.  range.  The  results  are  shown  in  Fig.  4.3. 
There  is  a  very  regular  decreasing  trend  in  the  matching  coefficient  as  range  increases.  The 
reasons  for  this  degradation  will  be  examined  in  the  next  subsection. 

Figures  4.4,  4.5,  and  4.6  show  transmission  loss  plots  and  Argand  diagrams  for 


Depth  (m)  Depth 
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Figure  4.3.  Averaged  extrapolation  results  V5.  range  for  Ex.  4. 

three  different  extrapolation  ranges.  The  matching  coefficients  and  ranges  are, 
respectively,  0.993  at  0.2  km,  0.948  at  3.0  km,  and  0.871  at  5.0  km.  As  before,  the  specific 
cases  shown  were  chosen  because  their  matching  coefficients  were  close  to  the  average 
matching  coefficient  for,  in  this  case,  the  given  range.  Notice  that  the  matching  coefficient 
for  0.2  km  is  only  slightly  below  that  for  the  unperturbed  case  (0.993  vs.  0.996).  Some 
other  interesting  characteristics  will  be  examined  in  the  next  subsection. 

4.5.3.  Discussion  of  Mismatched  Extrapolation 

We  can  analytically  examine  the  effects  of  the  water  column  sound  speed  profile 
mismatch  on  the  extrapolation  algorithm  by  going  back  to  our  normal  mode  representation 
for  the  modeled  field  (Eq.  2.5),  and  inserting  terms  for  the  perturbations: 


-il/2 


p(Zo) 


(c  +Ac  )  e 


(4.35) 


Multiplying  out,  ignoring  terms  that  multiply  perturbations  together,  and  expanding  the 


Depth  (m) 
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Held  extrapolation  at  3.0  km. 


Figure  4.6.  Ex.  4B:  Field  extrapolation  at  5.0  km. 
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inverse  square  root  of  the  perturbed  eigenvalue  in  a  first-order  Taylor  series  approximation 
results  in 


i  1 


87tr 


-1/2 


A4„5 


-3/2  \ 


e  e 


(4.36) 


This  expression  can  be  expanded,  once  again  ignoring  terms  that  multiply  perturbations 
together,  into  the  following  form: 


i 

1/2  1  p 

_ 1 

o 

Z  <'t'„(2o>AV„fe„)  +  AV„fe„)¥j2,))  . 


m  ^rn 


(4.37) 


In  Eq.  4.37,  we  take  note  of  the  first  summation  which  would  give  us  the 
expression  for  the  perfect  knowledge  modeled  field,  except  for  a  range-dependent  phase 
factor.  The  second  two  summations  have  perturbational  multipliers  (either  A\|/^(z)  or 
A^^).  For  Ex.  4B,  Fig.  4.2  and  Table  4.3  show  that  these  perturbations  are  very  small. 
Keeping  in  mind  that  this  expression  must  be  multiplied  by  the  Green’s  function  and 
integrated  over  depth,  we  also  take  note  that  for  short  ranges,  where  the  range-dependent 
phase  factors  will  have  little  effect,  the  matching  coefficient  is  very  close  to  that  of  the 
unperturbed  case.  This  implies  that  the  second  two  summations  in  Eq.  4.37  are  not 
exerting  a  meaningful  effect  on  the  outcome.  Since  the  effect  of  these  two  summations 
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will  not  vary  much  with  range,  we  will  neglect  their  contribution  to  the  modeled  field,  and 
write  the  averaged  extrapolated  field  as 


P'(^b’h')  = 


-1 1/2 


Snr. 


-1/2 


p(Zo) 


(4.38) 


We  are  assuming  that  for  this  scenario  the  obliquity  factor  did  not  appreciably  affect  the 
results. 

The  range-dependent  unit-amplitude  phase  factor  at  the  end  of  Eq.  4.38  is 
responsible  for  the  degradation  with  range  of  the  extrapolation  results  seen  in  Figs.  4.3-6. 
It  is  helpful  to  think  of  the  modeled  field  as  a  weighted  sum  of  the  mode  functions  as  in 
Eq.  1.11.  As  the  extrapolation  range  increases,  the  phase  of  this  term  moves  away  from 
zero,  causing  the  phase  of  the  complex  weight  for  each  mode  to  skew  from  its  “perfect 
knowledge”  position.  The  coherent  sum  of  the  weighted  mode  functions  is  then  changed. 
The  larger  the  extrapolation  range,  the  greater  the  phase  shift,  until  the  phase  wraps 
through  180°.  At  this  point  the  matching  coefficient  would  take  on  a  “random” 
characteristic,  depending  on  which  modes  had  phase  shifts  closer  to  zero  for  a  given  range. 

The  degradation  in  the  matching  coefficient  for  this  case  is  distinctly  different  from 
that  due  only  to  the  obliquity  factor  as  seen  in  Ch.  3.  The  obliquity  factors  are  amplitude 
factors  for  each  mode,  and  therefore  sum  in  phase  with  the  mode  they  are  related  to.  This 
effect  manifests  itself  in  the  transmission  loss  comparisons  in  Figs.  3.6-8  where  the 
transmission  loss  errors  tend  to  be  either  mostly  positive  or  mostly  negative,  depending  on 
whether  the  obliquity  factors  overweight  or  underweight  the  modes.  The  nature  of  the 
phase  errors  is  much  different  as  can  be  seen  in  Figs.  4.4-6.  Here  we  see  that  the 
transmission  loss  errors  are  more  evenly  distributed  about  0  dB. 

We  can  undertake  an  analysis  similar  to  Sec.  3.4  in  order  to  deterrmne  the  effect  of 
the  phase  multipliers  on  the  matching  coefficient.  We  now  write  the  fields  as 
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m 

P(Z)  =  J^W^\\fJz)e 


(4.39) 


Continuing  the  analysis  much  as  in  Eqs.  3.7-9,  we  arrive  at  the  relationship 


IFJ 


,2A2 


(4.40) 


Making  the  same  approximation  as  before,  namely  that  the  are  constant  with  m, 
allows  us  to  write  Eq.  4.40  as 


M 


2 


=  Y 


Ac  ’ 


(4.41) 


where  M  is  the  number  of  modes  that  contribute  significantly  to  the  field.  Note  that  for  all 
equal  to  zero  Y^^  is  equal  to  one,  which  would  be  the  expected  value  of  |isr|  for  this 

case. 

This  relation  implies  that,  if  one  has  estimates  of  the  changes  in  horizontal 
wavenumbers  due  to  perturbations  in  the  sound  speed  profile,  predictions  of  the  matching 
coefficient  can  be  obtained.  In  this  case,  it  is  useful  for  determining  how  changing  the 
variance  of  the  sound  speed  profile  perturbations  will  affect  the  matching  coefficient.  In 
Eqs.  4.31-33  we  saw  that,  for  the  case  of  additive  white  Gaussian  perturbations  to  the 
sound  speed  profile,  the  variance  of  the  noise  was  a  simple  multiplier  in  the  expression  for 
the  variance  of  the  perturbations  in  the  horizontal  wavenumbers. 
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Figure  4. 7.  Results  of  increased  variance  for  sound  speed  perturbations. 


If,  then,  for  the  same  realization  of  sound  speed  perturbations  as  seen  in  Fig.  4.1 
we  increase  the  variance  from  one  to  five,  each  horizontal  wavenumber  perturbation 
should  be  multiplied  by  Vs  s  2.236.  In  Eq.  4.41,  this  multiplier  could  be  factored  out  of 
the  eigenvalue  perturbation  and  would  become  a  multiplier  on  the  range,  meaning  that  the 
same  level  of  matching  coefficient  for  a  given  range  as  seen  in  Fig.  4.3  should  occur  at  that 
range  divided  by  approximately  2.236  for  the  increased  variance  case.  Figure  4.7  shows 
the  results  of  increasing  the  variance  in  the  sound  speed  perturbations  compared  to  the 
results  of  Fig.  4.3.  Note  that  for  the  first  case  the  matching  coefficient  is  approximately  0.6 
at  a  range  of  10  km,  and  at  the  same  level  at  approximately  4.4  km  for  the  increased 
variance  case  (10/2.236  =  4.472). 

4.6.  Mismatch  Effects:  Sediment  Parameters 

This  section  will  examine  the  effects  of  mis-estimation  of  some  of  the  sediment 
parameters  that  affect  acoustic  propagation,  namely  compressional  wave  (sound)  speed, 
compressional  wave  (sound)  attenuation,  and  density.  A  general  scenario  like  that  seen  in 
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Exs.  3F  and  4B,  i.e.,  a  sedimentary  layer  underlying  the  water  column,  will  be  used.  Note 
that  this  scenario  is  such  that  significant  propagation  of  acoustic  energy  in  the  sediment 
does  not  occur.  This  scenario  will  be  used  to  avoid  confounding  effects  in  the 
extrapolation  algorithm  due  to  truncated  arrays  since  we  are  assuming  that  the  arrays 
cannot  extend  into  the  bottom.  Effects  of  a  short  reference  array  in  deep  water  on  the  HAP 
algorithm  (see  Sec.  1.3)  were  examined  in  [Al-Kurd].  Bottom-propagating  acoustic  energy 
is  of  concern  in  some  underwater  scenarios,  and  its  effect  on  extrapolation  algorithms  is  an 
area  of  future  work. 

Accurate  estimation  of  bottom  parameters  that  affect  acoustic  propagation  is 
usually  a  more  difficult  task  than  that  of  determining  those  for  the  water  column.  The 
sound  speed  profile  in  the  water  column  is  the  most  important  property  there  and  is 
commonly  measured  with  CTD  (conductivity-temperature-depth)  instrumentation  [Wen  et 
al.].  Methods  for  determining  bottom  properties  range  from  coring  and  grab-sample 
methods  [Hamilton]  to  sophisticated  techniques  involving  mathematical  inversion  of 
acoustic  data  (see  articles  by  [Lindsay  and  Chapman]  and  [Dosso  et  al.},  and  references 
therein).  The  goal  of  this  section  is  to  attempt  to  gauge  the  relative  importance  of  the 
various  bottom  parameters  to  the  extrapolation  technique. 

4.6.1.  Sound  Speed  Mismatch  (Examples  4C  and  4D) 

For  Ex.  4C,  we  will  examine  the  effects  of  a  constant  offset  in  the  sound  speed 
estimate  in  the  sediment,  i.e.,  the  slope  of  the  sound  speed  profile  in  the  sediment  will 
remain  unchanged.  A  common  measure  for  describing  different  sediment  types  is  the 
sound  speed  ratio,  i.e.,  the  ratio  of  the  sound  speed  in  the  sediment  to  the  sound  speed  in 
the  water  at  the  sediment-water  boundary.  For  our  scenario  this  ratio  is 
1540/1480=  1.041 .  It  is  reported  in  [Hamilton]  that  sound  speed  ratios  can  vary  in  the 
range  from  0.976  to  1.201,  with  the  values  for  “clayey-silt”  compositions  lying  close  to 
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Table  4.4.  Eigenvalue  Comparison  for  Ex.  4C. 


Mode  #  (m) 

Eigenvalue  (^^) 

Change  in  Eigenvalue 

2 

1.267 +  i  1.171x10-5 

3.684x10-5  +  i  2.504x10-^ 

4 

1.262 +  i  1.121x10-5 

3.541x10-5  +  12.365x10-^ 

6 

1.259 +  i  1.166x10-5 

3.684x10-5  +  12.446x10-^ 

8 

1.254 +  i  1.696x10-5 

5.341x10-5  +  13.583x10-^ 

10 

1.248  +  12.530x10-5 

7.916x10-5  +  15.625x10-^ 

15 

1.228  +  18.751x10-5 

2.390x10-'^  +  i  3.449x10-5 

22 

1.197  +  15.809x10-'^ 

3.610x10-5-17.247x10-'^ 

29 

1.166 +  17.716x10-'^ 

3.573x10-5-16.682x10-'^ 

36 

1.136 +  i  1.110x10-5 

4.892x10-5-14.535x10-'^ 

43 

1.105 +  i  1.182x10-5 

4.022x10-5-12.841x10-'^ 

1 .0  (this  value  varies  in  the  vicinity  of  1 .0  for  clayey-silt  sediments  depending  on  the  type 
of  region  in  the  ocean  where  the  sediment  is  found,  e.g.,  continental  terrace,  abyssal  hill, 
or  abyssal  plain). 

Example  4C  will  use  a  constant  offset  of  -10  m/s  in  the  sediment  sound  speed 
resulting  in  a  sound  speed  of  1530  m/s  at  the  sediment/water  interface.  All  other 
parameters  are  the  same  as  Ex.  3F  (including  the  slope  of  the  sound  speed  profile  in  the 
sediment).  Note  that  we  are  assuming  perfect  knowledge  of  the  sound  speed  profile  in  the 
water  column.  Table  4.4  shows  a  subset  of  the  computed  eigenvalues  and  their  change 
from  the  unperturbed  case.  Notice  the  increasing  trend  in  the  perturbations  to  the 
eigenvalues,  and  that  all  the  perturbations  to  the  real  parts  of  the  eigenvalues  are  positive. 
We  can  use  the  results  derived  in  Sec.  4.3  to  explain  this  phenomenon. 

Even  though  Eq.  4.23  was  derived  while  considering  sound  speed  perturbations  in 
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the  water  column,  it  is  also  valid  for  sound  speed  perturbations  in  the  sediment.  Since  we 
are  assuming  no  perturbations  in  the  water  column,  the  integral  now  spans  h  (instead  of 
zero)  to  oo  and  Ac(z')  becomes  a  constant  with  depth  that  can  be  taken  outside  the 
integral.  Using  the  same  reasoning  found  in  Sec.  4.5.1,  we  can  see  that  modes  whose  mode 
functions  have  more  “energy”  in  the  sediment  will  be  affected  to  a  greater  extent  by  the 
perturbations  to  the  sound  speed  estimate  in  the  sediment.  Also,  since  the  perturbation  to 
the  sound  speed  is  a  negative  constant,  it  comes  out  of  the  integral  and  multiplies  the  -1 
already  in  front  of  the  integral  to  create  a  positive  multiplier.  The  rest  of  the  integrand  is 
then  positive,  so  the  perturbations  to  the  real  parts  of  the  eigenvalues  are  positive  in  this 
case  (the  background  wavenumber  and  sound  speed  can  actually  be  thought  of  as  complex 
numbers  due  to  the  attenuation  in  the  sediment,  as  will  be  discussed  in  the  next 
subsection). 

Figure  4.8  shows  a  subset  of  the  perturbed  (black  trace)  vs.  unperturbed  (gray 
trace)  mode  functions.  Notice  that  the  lower  order  (propagating,  in  this  case)  modes  have 
no  noticeable  perturbations.  Again  Eq.  4.26  gives  us  a  clue  as  to  why  this  is  the  case.  The 
integrals  now  need  only  cover  the  sediment  region  since  the  sound  speed  profile  in  the 
water  column  is  unperturbed.  Since  tbe  lower  order  modes  have  little  extension  into  the 
sediment,  their  perturbations  are  small. 

Figure  4.9  shows  some  extrapolation  results  for  this  mis-estimation  of  the 
sediment  sound  speed  profile.  It  is  a  plot  of  matching  coefficient  vs.  extrapolation  range, 
plotted  on  the  same  scale  as  Fig.  4.3  for  comparison.  We  can  use  Eq.  4.41  to  explain  the 
much  reduced  falloff  compared  to  the  perturbed  water  column  sound  speed  profile  case 
(Ex.  4B).  Note  that  the  perturbations  to  the  eigenvalues  (for  the  non-stripped  modes)  are 
roughly  of  the  same  order  of  magnitude  as  those  for  the  Ex.  4B,  even  though  the  sediment 
mismatch  is  10  m/s,  whereas  the  random  perturbations  to  the  water  column  sound  speed 
profile  had  a  variance  of  only  1  m/s.  Also  note  that  the  eigenvalue  perturbations  for  the 
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Figure  4.8.  Comparison  of  “true"  and  perturbed  mode  functions  for  Ex.  4C 
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Figure  4.9.  Ex.  4C:  Constant  ojfset  in  sediment  sound  speed  profile. 

water  column  fluctuations  were  not  consistent  in  sign.  In  Eq.  4.41,  we  see  that  is  a 
function  of  the  summation  of  the  exponentiated  eigenvalue  perturbations.  When  the 
are  all  the  same  sign,  they  will  sum  in  phase  (assuming  rA^^  «  n  for  all  m,  which  is  the 
case  for  ranges  considered  here)  as  opposed  to  the  case  where  they  are  of  differing  sign; 
hence  the  results  seen  in  Fig.  4.9. 

For  a  more  direct  comparison,  we  will  perturb  the  sound  speed  in  the  sediment 
with  the  same  random  fluctuations  that  were  added  to  the  water  column  sound  speed 
profile  (Ex.  4D).  Table  4.5  shows  the  eigenvalues  for  this  perturbed  system  and  the 
changes  from  the  unperturbed  state.  In  comparison  to  Table  4.3  (water  column  sound 
speed  fluctuation  case),  we  see  that  the  perturbations  to  the  real  parts  of  the  eigenvalues 
are  roughly  an  order  of  magnitude  less  and  of  consistent  sign.  In  comparison  with  Table 
4.4  (constant  offset  in  sediment  sound  speed  profile),  we  see  that  the  perturbations  are  also 
roughly  an  order  of  magnitude  less.  We  would  expect,  then,  that  the  degradation  of  the 
extrapolated  fields  as  expressed  in  the  matching  coefficient  would  be  much  less  for  Ex.  4D 
than  for  Exs.  4B  or  4C.  This  is  indeed  the  case  as  seen  in  Fig.  4.10  which  shows  the 
matching  coefficient  vs.  extrapolation  range  for  Ex.  4D.  As  seen  here,  the  matching 


Table  4.5.  Eigenvalue  Comparison  for  Ex.  4D. 


Mode  #  (m) 

Eigenvalue 

Change  in  Eigenvalue 

2 

1.267  +  19.115x10'^ 

-1.311x10'^ -18.599x10'* 

4 

1.262  +  18.760x10'^ 

-1.192x10'^ -18.535x10'* 

6 

1.259  +  19.125x10'^ 

-1.311x10'^ -19.294x10'* 

8 

1.254 +  i  1.323x10'^ 

-2.027x10'^  - 1  1.456x10''^ 

10 

1.248 +  i  1.942x10'^ 

-2.981x10-^-12.520x10''^ 

15 

1.228  +  i  5.084x10'^ 

-1.204x10-^-12.187x10'^ 

22 

1.191  +18.027x10''^ 

-1.740x10-^-15.028x10''^ 

29 

1.158  +  1  1.035x10'^ 

-3.960x10-^-14.044x10''^ 

36 

1.126  +  19.483x10'^ 

-4.287x10-3-16.150x10''^ 

43 

1.095 +  i  1.530x10'^ 

-5.551x10-3-16.405x10-3 

Figure  4.10.  Ex.  4D:  Random  fluctuations  in  sediment  sound  speed  profile. 

coefficient  hardly  deviates  at  all  from  its  ideal  value  of  1 .0. 

The  general  conclusion  we  can  draw  from  the  above  results  is  that  for 
environmental  scenarios  that  do  not  support  significant  propagation  of  sound  through  the 
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sediment,  the  most  important  part  of  estimating  the  sound  speed  profile  in  the  sediment  is 
in  the  upper  most  regions  of  the  sediment  near  the  sediment/water  interface  where 
boundary  conditions  are  established.  While  intuitively  satisfying,  we  can  also  see  why  this 
is  the  case  due  to  our  analysis  of  the  perturbations  to  the  eigenvalues.  The  relationships 
governing  these  perturbations  involve  integrals  over  the  mode  functions  where  the  sound 
speed  perturbations  exist.  If  the  mode  functions  are  small  where  the  sound  speed 
perturbations  exist,  then  the  perturbations  to  the  horizontal  wavenumbers  (and  mode 
functions)  will  be  small. 

4.6.2.  Attenuation  Mismatch  (Example  4E) 

The  second  sediment  parameter  whose  effects  on  the  extrapolation  we  are 
interested  in  gauging  is  attenuation.  As  was  seen  in  Ch.  3,  imaginary  components  of  the 
horizontal  wavenumbers  manifest  themselves  in  mode-dependent  attenuations.  The  usual 
analysis  of  attenuation  [Jensen  et  al,  Boyles]  begins  by  assuming  a  small  imaginary 
perturbation  to  the  real-valued  sound  speed  profile,  i.e., 

c(z)  =  c^(z)-ic.(z).  (4.42) 

This  perturbation  is  negative  because  when  the  attenuation  rule  in  general  is  written  as 
^(ikA  a.r)  ^  ^  wavenumber  and  a  is  in  units  of  nepers/m,  we  can  set 

this  equal  to 


i(Ox 

C 

e 


iCOA 

e  ^ 


Cr  + 


>  2 
■  +  C; 


(4.43) 


assuming  that  the  squared  c^.’s  are  negligible.  This  allows  us  to  make  the  approximation 
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(4.44) 


meaning  that  a  positive  value  of  c  •  leads  to  a  positive  value  of  a . 

We  can  now  undertake  a  perturbational  analysis  much  as  in  Sec.  4.3  in  order  to 
determine  what  the  perturbations  to  the  horizontal  wavenumbers  will  be  for  a  perturbation 
in  attenuation  [Rajan  et  al].  Using  Eq.  4.42  as  our  perturbation  and  starting  with  a 
relationship  much  like  Eq.  4.16,  we  eventually  arrive  at 


(1) 

m 


(0), 


2C.(z') 


(0) 
m  0 


^o(^) 


(4.45) 


Note  that  this  is  positive  and  wholly  imaginary-valued.  Substituting  in  for  c-  via  Eq.  4.44, 
we  see  that 


0 

where  a  is  in  units  of  nepers/m.  The  model  we  have  been  using  expressed  its  attenuation 
in  units  of  dB/wavelength.  Using  the  approximation  that  a  «  8.686aX  [Jensen  et  a/.], 
we  can  write  the  perturbations  to  the  horizontal  wavenumbers  as 

oo 

2,t  (8.686) 

where  is  now  in  units  of  dB/wavelength.  Since  this  analysis  started  from  the 

assumption  of  no  attenuation  (real-valued  sound  speed  profile),  determining  the  change  in 
eigenvalues  where  the  starting  scenario  had  attenuation  would  require  computations  of  the 
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eigenvalue  changes  via  Eq.  4.47  for  both  scenarios  and  a  differencing  operation. 

The  perturbations  to  the  mode  functions  are  wholly  imaginary-valued  as  well. 
These  perturbations  would  be  very  small,  and  KRAKEN  ignores  them  for  the  modes  it 
calculates,  so  we  will  neglect  them  for  this  analysis. 

We  will  again  use  the  environmental  scenario  of  Ex.  3F  as  our  ground  truth  case 
(see  Fig.  3.10),  and  now  assume  a  mis-estimation  of  the  sediment  attenuation  value.  As 
reported  in  [Hamilton],  attenuation  values  vary  roughly  over  a  range  from  very  near  0.0  to 
1.275  dB/A,.  Also,  it  is  mentioned  in  [Rajan  et  a/.]  that  it  is  common  practice  to  assume  a 
constant  attenuation  value  with  depth,  as  was  done  for  our  model.  The  perturbation  here 
(Ex.  4E)  will  consist  of  mis-estimating  the  sediment  attenuation  as  0.2  dB/A,  vs.  the 
ground  truth  value  of  0. 1  dB/A . 

Table  4.6  shows  the  eigenvalues  and  the  changes  from  the  ground  truth  case  for 
Ex.  4E.  Note  that  while  there  are  some  perturbations  to  the  real  parts  of  the  eigenvalues, 
they  are  extremely  small  compared  to  the  actual  eigenvalues,  and  hence  we  will  neglect 
them.  The  perturbations  to  the  imaginary  parts,  however,  are  comparable  in  order  of 
magnitude  to  the  base  values.  Noting  that  the  imaginary  parts  of  the  eigenvalues  are 
responsible  for  attenuations  imparted  to  each  mode,  we  can  compare  the  attenuations  over 
a  20  km  range  due  to  the  ground  truth  and  perturbed  eigenvalues  (see  Fig.  3.12)  as  in  Fig. 
4.11.  The  perturbed  results  are  in  gray.  As  we  would  expect,  the  higher  attenuation  value 
has  resulted  in  greater  attenuation  being  imparted  to  each  mode. 

Extrapolation  results  (matching  coefficient  vs.  range)  for  this  case  are  shown  in 
Fig.  4.12.  The  degradation  in  the  extrapolated  field  as  expressed  by  the  matching 
coefficient  is  quite  small  for  this  case. 

In  general,  since  the  perturbations  to  the  imaginary  parts  of  the  eigenvalues 
manifest  themselves  in  modal-dependent  amplitude  factors,  we  can  use  an  analysis  similar 
to  that  for  the  obliquity  factors  to  derive  an  expression  that  predicts  the  degradation  in 
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Table  4.6.  Eigenvalue  Comparison  for  Ex.  4E. 


Mode  #  (m) 

Eigenvalue 

Change  in  Eigenvalue 

2 

1.267 +  i  1.839x10'^ 

0  +  19.192x10-^ 

4 

1.262 +  i  1.768x10-^ 

0  +  18.837x10-^ 

6 

1.259 +  i  1.843x10*^ 

0  +  19.209x10-^ 

8 

1.254  +  12.674x10-^ 

-1.190x10-'^  +  !  1.336x10-2 

10 

1.248  +  13.933x10-^ 

-1.200x10-'^  +  !  1.965x10-2 

15 

1.228 +  i  1.060x10-4 

-4.770x10-^  +  i  5.293x10-2 

22 

1.193  +  12.609x10-2 

-8.703x10-^  +  1  1.304x10-2 

29 

1.162  +  12.882x10-2 

-9.298x10-^  +  1  1.442x10-2 

36 

1.131+13.125x10-2 

-9.537x10-^  +  1  1.562x10-2 

43 

1.101  +12.935x10-2 

-9.179x10-^  +  1  1.469x10-2 

Figure  4.11.  Ex.  4E:  Attenuation  due  to  imaginary  part  of  eigenvalues. 


matching  coefficient  for  given  perturbations  in  the  imaginary  parts  of  the  eigenvalues.  In 
this  case  we  arrive  at  the  relationship 
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Figure  4.12.  Ex.  4E:  Constant  offset  in  sediment  attenuation. 
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(4.48) 


4.6.3.  Density  Mismatch  (Example  4F) 

It  is  shown  in  [Raj an  et  al]  that  similar  results  can  be  obtained  for  perturbations  in 
the  density  profile  as  were  obtained  for  perturbations  to  the  sound  speed  and  attenuation. 
Specifically,  we  can  define  an  intermediate  function 


m{z)  = 


p'iz) 

2 


from  which  the  perturbations  to  the  eigenvalues  can  be  written  as 


(4.49) 


(1) 


1 


0 


J(Po(z'))' 


^Mz')dz' . 


(4.50) 


The  perturbations  to  the  mode  functions  follow  in  a  similar  manner. 
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Table  4.7.  Eigenvalue  Comparison  for  Ex.  4E. 


Mode  #  (m) 

Eigenvalue 

Change  in  Eigenvalue 

2 

1.267 +  i  1.029x10'^ 

1.162x10''^ +  i  1.092x10'^ 

4 

1.263  +  19.065x10'^ 

9.620x10'^ +  12.196x10''^ 

6 

1.259  +  18.875x10'^ 

8.798x10'^  -  i  3.432x10'^ 

8 

1.254 +  i  1.219x10'^ 

1.098x10'^ -i  1.191x10'^ 

10 

1.249 +  i  1.702x10'^ 

1.305x10'^ -12.652x10'^ 

15 

1.228  +  14.092x10'^ 

9.859x10'^  -  i  1.210x10'^ 

22 

;  1.193  +  1  1.448x10'^ 

-6.414x10'^  +  i  1.419x10'^ 

29 

1.162  +  1  1.561x10'^ 

8.893x10'^ +  i  1.209x10''^ 

36 

1.131 +i  1.675x10'^ 

-7.260x10'^ +  i  1.117x10'^ 

43 

1.101  +i  1.519x10'^ 

1.695x10'^ +  15.254x10'^ 

Values  of  density  found  in  sediments  vary  from  lows  of  roughly  1.2  g/cm^ 
[Hamilton]  to  2.7  g/cm^  [Jensen  et  al.].  For  Ex.  4F  we  will  perturb  our  ground  truth 
estimate  by  changing  the  density  to  2.0  g/cm^  vs.  the  true  value  of  1.5  glen?.  Table  4.7 
shows  a  subset  of  the  perturbed  eigenvalues  and  their  change  from  the  ground  truth  case. 
The  perturbed  mode  functions  showed  essentially  no  change  from  the  ground  truth,  so 
they  are  not  shown. 

Figure  4.13  shows  the  matching  coefficient  vs.  range  results  for  Ex.  4F.  Since  the 
density  perturbations  affect  mainly  the  real  part  of  the  eigenvalues  (like  sound  speed 
perturbations),  we  can  use  the  expression  in  Eq.  4.41  to  predict  the  degradation  in  the 
matching  coefficient  due  to  the  perturbations  in  the  density,  i.e.,  Y^p  =  Y^^. 

4.6.4.  Discussion  of  Sediment  Perturbation  Results 

The  goals  of  this  section  were  twofold.  One  was  to  reinforce  the  finding  in 
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Figure  4.13.  Ex.  4F:  Constant  offset  in  sediment  density. 

[Rouseff  et  al]  that  for  environments  where  sediment-borne  energy  is  negligible,  the 
parameters  in  the  bulk  of  the  sediment  that  affect  propagation  are  significantly  less 
important  relative  to  other  factors  such  as  the  sound  speed  in  the  water  column  or  the 
bathymetry  (the  bathymetry  issue  was  dealt  with  in  [Rouseff  et  al.]).  The  examples  given 
here  showed  that  for  relatively  large  perturbations  to  the  sediment  parameters,  the 
degradation  in  matching  coefficient  was  significantly  less  than  that  caused  by 
perturbations  to  the  water  column  sound  speed  profile.  These  differences  were  explained 
by  appealing  to  the  perturbational  theory  results  that  showed  how  modes  were  affected  by 
the  parameter  perturbations  as  a  function  of  how  far  their  mode  functions  extended  into 
the  sediment.  In  this  case,  where  the  mode  functions  that  contributed  significantly  to  the 
fields  did  not  extend  very  far  into  the  sediment,  the  perturbations  were  relatively  small  and 
the  extrapolated  fields  exhibited  high  matching  coefficients  relative  to  the  ground  truth 
fields. 

The  second  goal  was  to  give  an  approximate  measure  by  which  one  could,  given  a 
certain  perturbed  scenario,  determine  how  the  matching  coefficient  would  be  affected. 
This  was  accomplished  through  the  T  expressions.  In  the  preceding  text,  individual  T ’s 
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were  derived  depending  on  whether  the  perturbations  to  the  modes  were  phase  or 
amplitude  factors.  We  note  in  summary  that  a  combined  T  can  be  defined  by  assuming  the 
extrapolated  field  is  compose  of  modes  with  a  complex  weighting  factor,  i.e., 

(4.51) 

m 

This  results  in  the  more  general  relation 


Y  = 


M 


m  =  1 


M 


m=  1 


(4.52) 


4.7.  Summary  of  Chapter  4 

Chapter  4  dealt  with  the  effects  on  the  extrapolation  algorithm  of  mis-estimating 
some  of  the  environmental  parameters  that  are  used  as  input  to  the  modeled  field 
processing.  Expressions  were  derived  which  showed  how  perturbations  to  environmental 
parameters  affect  the  mode  functions  and  horizontal  wavenumbers  of  the  modeled  field. 
These  expressions  were  then  interpreted  with  the  aid  of  several  example  scenarios.  It  was 
seen  for  scenarios  where  a  sediment  layer  strips  off  the  higher  order,  more  steeply 
propagating  modes,  that  the  extrapolation  algorithm  degrades  more  rapidly  with 
extrapolation  range  due  to  perturbations  in  the  water  column  sound  speed  profile  than  due 
to  relatively  larger  perturbations  to  sediment  parameters  (sound  speed,  attenuation,  and 
density).  Expressions  were  shown  which  allow  estimation  of  the  matching  coefficient 
given  estimates  of  the  perturbations  to  the  horizontal  wavenumbers.  Chapter  5  will 
explore  numerical  integration  techniques  for  the  extrapolation  algorithm. 


Chapter  5:  Integral  Estimation  for  Acoustic  Field  Extrapolation 


5.1.  Introduction  to  Chapter  5 

In  Ch.  2,  the  basic  theory  behind  the  acoustic  field  extrapolation  algorithm  was 
explored.  There  we  saw  that  calculating  the  integral  of  the  product  of  either  the  modeled 
field  or  its  derivative  and  the  measured  calibration  data  (i.e.,  Green’s  function  data)  is 
central  to  the  algorithm.  Practicality  requires  that  the  calibration  data  be  collected  for 
discrete  transmitter  and  receiver  locations  along  the  reference  arrays,  meaning  that  a 
numerical  integration,  or  “quadrature,”  technique  must  be  used  to  estimate  the  integral. 

This  chapter  will  explore  issues  related  to  estimating  the  integrals  as  in  Eq.  2.4.  It 
will  begin  with  a  discussion  of  trapezoidal  integration,  and  how  the  extrapolation 
algorithm  degrades  as  the  density  of  transmitters  along  the  reference  array  decreases  when 
this  numerical  integration  method  is  used.  It  will  suggest  methods  by  which  we  can 
estimate  the  spacing  of  elements  required  along  the  reference  arrays  in  order  to  extrapolate 
fields  reliably.  Another  numerical  integration  technique,  called  an  inner  product 
quadrature  formula,  will  be  examined  which  uses  the  modeled  mode  functions  as  an 
integration  basis.  This  integration  method  will  show  the  relationship  between  the 
Huygens’  principle-related  extrapolation  algorithm  and  a  modal  decomposition-based 
method. 

For  this  material  we  will  assume  that  the  reference  transmit  array  spans  the  water 
column  and  the  transmitters  are  spaced  equally  along  the  array.  Note  that  for  this 
application  our  degrees  of  freedom  are  restricted  in  that  the  abscissa  points  where  the 
function  to  be  integrated  is  evaluated  are  defined  by  the  locations  in  depth  of  the  reference 
transmitters.  This  means  that  integration  methods  that  have  flexibility  in  choosing  the 
abscissa  locations,  thereby  adding  more  degrees  of  freedom,  cannot  be  used  (e.g.,  certain 
types  of  Gaussian  quadratures  [Press  et  al,  Davis  and  Rabinowitz].) 
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5.2.  Trapezoidal  Integration 
5.2.1.  Definition 

A  common  and  easily  implementable  method  of  numerical  integration  for  evenly 
spaced  sample  points  is  trapezoidal  integration  [Press  et  al].  Trapezoidal  integration  is 
performed  by  assuming  straight  lines  between  the  ordinates  of  the  function  to  be 
integrated,  and  summing  the  areas  of  the  resulting  trapezoids  as  shown  in  Fig.  5.1.  For  this 
example  we  can  say  that 


^y{x)dx~  ^  areas  of  trapezoids . 
0 


The  area  of  any  given  trapezoid  is 


A 


n 


h - r - 


(5.1) 


(5.2) 


By  expanding  the  expressions  for  all  the  trapezoids  and  collecting  like  terms,  we  find  that 
the  integral  of  the  function  can  be  approximated  by 


Figure  5.1.  Example  of  trapezoidal  integration. 
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I  y  (x)  dx~h 


'^  +  XD  +  >(2)  +  J(3)  +  y(4)  +  ^' 


(5.3) 


The  general  trapezoidal  integration  rule  is  given  by 


J  >•  (x)  dx  -  h\^  +  y(x,)  +  yix^)  +  . . .  +  yix^  _  2)  +  y(x^  _  j)  +  J.  (5.4) 

*0 

Note  that  except  for  the  end  conditions,  this  equates  to  simply  summing  the  samples  of  the 
function  to  be  integrated  and  multiplying  by  the  interval  between  samples. 

5.2.2.  Fourier  Theory  Considerations 

We  define  the  Fourier  transform  of  a  continuous-time  function  x{t)  as 


X(f)  =  jx(t)e~-’^''^‘dt. 


where  f  is  frequency.  This  leads  to  the  expression 


(5.5) 


X(0)  =  J  x{t)dl.  (5.6) 

_CXD 

In  other  words,  the  integral  of  a  function  over  all  time  is  equal  to  the  “DC”  term  of  its 
Fourier  transform. 

If  we  assume  that  a  given  x(r)  is  bandlimited,  i.e.,  its  Fourier  transform  is 
identically  zero  for  all  /|  >  fi,  then  using  a  Fourier  series  approach  [Marks]  we  can  arrive 
at  the  relationship 


Ill 


OO  oo 

j  x{t)dt  =T  Y,  x(nT),  (5.7) 

-oo  n  =  -oo 

where  T  is  the  sampling  interval  and  l/T>2B  in  order  to  satisfy  the  Nyquist  sampling 
criterion.  Note  the  similarity  to  the  trapezoidal  integration  rule,  Eq.  4.4.  Since  a  perfectly 
bandlimited  function  must  exist  for  all  time  -<»  <  r  <  oo ,  the  end  points  do  not  affect  the 
result.  We  see  then  that  a  properly  sampled  bandlimited  function  can  be  integrated  exactly 
using  trapezoidal  integration. 

Signals  sampled  below  the  Nyquist  frequency  can  also  be  integrated  exactly. 
Figure  5.2  illustrates  why.  The  top  picture  in  Fig.  5.2  shows  a  properly  sampled 
bandlimited  signal  with  l/Tj  >  25  so  that  there  is  no  aliasing  of  the  replicated  spectra. 
The  bottom  picture  shows  the  results  of  sampling  with  5  <  1/72  <25. 
show  where  the  baseband  spectrum  is  contaminated  by  the  overlapping  replicas  due  to 


25 
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Figure  5.2.  Frequency  domain:  aliasing  of  sampled  signals. 
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aliasing.  Note  that  the  DC  value,  X(0),  is  not  contaminated.  Therefore,  exact  integration 
results  can  be  obtained  for  bandlimited  signals  by  summing  the  samples,  assuming  that  the 
sampling  interval  is  such  that  B  <\/T  <1B . 

Truly  bandlimited  signals  (which  must  exist  for  all  time)  are  mostly  mathematical 
niceties.  Many  real-world  signals,  however,  are  “essentially”  bandlimited  in  that  their 
energy  content  above  a  given  frequency  limit  is  very  small.  Since  these  signals  are  not 
truly  bandlimited,  Eq.  5.7  is  not  an  equality.  As  the  sampling  interval  decreases  (sampling 
rate  increases),  however,  and  the  essential  bandwidth  of  the  signal  is  accounted  for  in  the 
sampling,  Eq.  5.7  becomes  a  better  approximation.  Using  the  above  ideas  we  can 
determine  the  sampling  interval  at  which  a  trapezoidal  integration  will  start  to  fail  to 
provide  a  close  estimate  of  the  true  integral. 

5.2.3.  Trapezoidal  Integration  Example 

For  this  example  we  will  use  the  same  environmental  scenario  as  Ex.  3F 
(attenuating  sediment  layer  underlying  the  water  column)  so  that  obliquity  factor  effects 
are  negligible.  In  Section  3.4.2  we  saw  that  using  a  reference  transmit  array  whose 
transmitters  were  spaced  1  m  apart  resulted  in  a  matching  coefficient  of  approximately 
0.996  for  an  extrapolation  from  =  (5  km,  30  m).  Figure  5.3  shows  the 

magnitude-squared  of  the  discrete  Fourier  transforms  (DFT’s,  calculated  by  zero  padding 
out  to  1024  and  using  an  FFT  routine,  see  [Oppenheim  and  Schafer])  of  a)  the  modeled 
field  p(5,  zj,  b)  a  representative  Green’s  function,  specifically  that  for  a  receive  depth  of 
74  m  G(25,  74  ;5,  zj ,  and  c)  the  product  of  the  modeled  field  and  Green’s  function.  It  is 
important  to  remember  that  in  the  acoustic  field  extrapolation  algorithm  an  integral  must 
be  approximated  for  each  point  along  the  receive  array  where  an  extrapolated  field  value  is 
desired.  The  Green’s  function  in  this  case  is  the  result  of  transmitting  from  each  of  the 
transmitters  along  the  reference  array  and  receiving  at  depth  74  m  at  the  receive  array. 
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hence  the  depth  variable  .  Using  the  principle  of  acoustic  reciprocity  (see  Section  2.5.8), 
this  can  be  seen  as  equivalent  to  transmitting  from  the  receive  location  and  receiving  along 
the  reference  transmit  array. 

Figure  5.3a  shows  the  DFT  of  the  modeled  field.  The  horizontal  axis,  or  what  we 
would  think  of  as  the  frequency  axis  (units  Hz,  or  1/s)  for  the  DFT  of  a  time-domain 
signal,  is  in  normalized  units  of  1/m;  0.5  units  on  this  axis  refers  to  one-half  the  spatial 
sampling  “rate.”  The  gray  vertical  markers  show  an  “eyeball”  estimate  of  the  essential 
“bandwidth”  of  the  sampled  field.  Figure  5.3b  shows  the  DFT  of  the  Green’s  function  for 
the  stated  receive  depth.  Note  that  the  essential  bandwidth  estimate  appears  very  close  to 
that  of  the  modeled  field.  Keeping  in  mind  the  reciprocity  relationship  for  the  Green’s 
function  and  interpretation  for  the  fields  based  on  the  sum  of  weighted  mode  functions,  we 
find  this  result  intuitively  satisfying  since  both  functions  are  weighted  sums  of  the  same 
mode  functions. 

Figure  5.3c  shows  the  DFT  of  the  product  of  the  modeled  field  and  the  Green’s 
function.  The  essential  bandwidth  here  is  roughly  twice  that  seen  in  Figs.  5.3a  and  b. 
Again,  this  result  is  intuitively  satisfying  since  multiplying  functions  in  the  spatially 
sampled  domain  (equivalent  to  the  temporally  sampled  domain  for  functions  of  time) 
translates  to  convolving  their  representations  in  the  spatial  “frequency”  domain  (often 
referred  to  as  the  modulation  property;  see  [Oppenheim  et  al.]). 

We  can  now  predict  the  effect  that  increasing  the  spacing  of  the  elements  on  the 
reference  array  (decreasing  the  spatial  sampling  “rate”)  has  on  our  extrapolation 
algorithm.  In  Fig.  5.3c  we  see  that  the  spectral  content  is  not  zero  outside  the  essential 
bandwidth  region,  but  rather  gradually  falls  off.  From  the  results  of  Sec.  5.2,  we  can 
expect  that  as  the  spacing  is  increased  from  1  m,  the  extrapolation  results  will  degrade  as 
this  energy  is  aliased  into  X  (0) .  When  the  outside  edges  of  the  “essential”  bandwidth 
region  (marked  by  the  bandwidth  delimiters  in  Fig.  5.3c)  begin  to  be  aliased  into  X  (0) , 
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Figure  5.4.  Degradation  of  matching  coefficient  with  increased  array  spacing. 

we  should  expect  a  more  drastic  degradation  of  the  extrapolation  results.  From  Fig.  5.3c, 
we  see  that  the  approximate  bandwidth  delimiter  is  placed  at  roughly  0.11  normalized 
units.  The  factor  that  would  cause  this  portion  to  begin  to  be  aliased  into  X(0)  is 
1  /0. 1 1  =  9.09 .  Therefore,  we  can  expect  more  serious  degradation  to  occur  near  spacings 
of  approximately  9.09  x  Im  =  9.09m. 

Figure  5.4  shows  the  results  of  varying  the  spacing  of  elements  along  the  reference 
transmit  array  for  the  previously  stated  scenario.  As  in  previous  simulations,  the  matching 
coefficient  results  have  been  averaged  over  a  grid  of  source  locations.  This  plot  shows  the 
anticipated  mild  degradation  of  the  results  to  a  certain  spacing,  followed  by  a  sharp 
downturn.  Note  that  the  predicted  location  of  the  sharp  downturn,  -9.09  m,  does  occur  in 
the  “knee”  of  the  curve. 

This  analysis,  then,  suggests  a  method  by  which  we  could  estimate  the  spacing 
needed  for  a  reference  array  in  a  given  environment.  Assuming  that  reasonably  reliable 
modeled  fields  can  be  produced  (which  would  be  necessary  for  the  success  of  the 
extrapolation  algorithm  in  any  case),  we  could  produce  a  set  of  modeled  fields  for  a 
representative  grid  of  hypothetical  source  locations.  Noting  that  spacing  of  the  modeled 
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field  points  in  depth  is  not  restricted,  we  could  then  make  an  estimate  of  the  essential 
bandwidth  of  the  modeled  fields  by  discrete  Fourier  transforming  in  depth.  The  grid  of 
source  locations  should  produce  a  good  idea  of  the  maximum  essential  bandwidth.  If  no  a 
priori  knowledge  of  the  essential  bandwidth  of  the  Green’s  function  is  available,  merely 
doubling  the  bandwidth  obtained  for  the  modeled  fields  can  serve  as  a  first-order 
approximation  of  the  product  bandwidth.  Otherwise,  information  about  the  bandwidth  of 
the  Green’s  function  can  be  incorporated  into  the  estimate.  Using  an  analysis  similar  to 
that  stated  above,  it  can  be  estimated  at  what  spacing  the  essential  bandwidth  portion  of 
the  product  Fourier  transform  will  begin  interfering  with  the  X(0)  term  and  drastically 
reducing  the  matching  coefficient.  As  can  be  seen  from  Fig.  5.4,  decreasing  the  spacing 
(adding  more  transmitters)  will  gradually  aid  the  performance  of  the  extrapolation 
algorithm. 

5.3.  Inner  Product  Quadrature  Formulas  for  Numerical  Integration 

Some  types  of  numerical  integration  methods  involve  fitting  functions  with 
integrals  that  can  be  calculated  analytically  to  the  sampled  values  of  the  function  to  be 
integrated.  Many  of  the  methods  use  orthogonal  polynomials,  such  as  Legendre  or 
Chebyshev,  for  the  basis  functions  [Press  et  al.].  Noting  for  now  that  we  can  produce  a 
“customized”  basis  set  of  orthonormal  functions  via  a  numerical  normal  mode 
decomposition,  we  will  now  outline  the  basic  theory  behind  an  integration  method  through 
which  we  can  take  advantage  of  this  fact.  In  the  following  section  some  results  from 
[Gribble]  and  [McGrath]  are  summarized. 

5.3.1.  Theory  for  Inner  Product  Quadrature  Formulas 

The  trapezoidal  rule  is  one  example  of  a  so-called  Regular  Quadrature  Formula 


(R.Q.F),  i.e., 
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f(z)dz = /(/)  ~  e  =  X 

j  =  0 

where  /?  is  a  real  bounded  interval  and  the  a. ’s  are  coefficients  that  define  the  particular 
R.Q.F.  For  integrands  that  are  inherently  products  of  functions  (as  is  our  field 
extrapolation  algorithm),  i.e., 

j wix)f(x)g(x)dx  =  I(f;g) ,  (5.9) 

R 

where  w  (jc)  is  a  nonnegative  weight  function,  there  exist  what  are  called  Inner  Product 
Quadrature  Formulas  (I.P.Q.F.’s). 

In  general,  these  I.P.Q.F.’s  take  the  form 

m  n 

Kf\g)  ~  Q  ” V;^)  =  E  S  (5-10) 

i  =  0j  =  0 

where  the  y.’s  and  Zj's  are  the  nodes  at  which  f  and  g  respectively  are  evaluated,  and  the 
a^j’s  are  determined  by  requiring  Eq.  5.10  to  be  exact  when  f  and  g  are  functions  that  are 

linear  combinations  of  sets  of  linearly  independent  basis  functions,  i.e., 

-  {5o>  •••>  ^5^  •  (5.11) 

Eq.  5.10  can  be  equivalently  expressed  in  matrix  form  as 

(5.12) 


Note  that  f  and  g  are  treated  individually  in  that  they  can  be  evaluated  on  separate  sets  of 
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nodes  (i.e.,  {Xg, x^}  and  {jg, not  coincident)  or  have  different  linear  bases 
if  need  be. 

To  solve  for  the  matrix  A ,  start  by  defining  the  matrices 


F  =  {/,(yp} 
G  =  {g,(z,.)}' 


1, 7  =  0 


i,  7  =  0 


(5.13) 


so  that  each  row  of  the  matrices  contains  the  samples  of  a  given  basis  function  evaluated  at 
the  nodes  assigned  to  that  basis  as  defined  in  Eq.  5.10.  According  to  [Gribble]  and 
[McGrath],  the  algorithm  proceeds  by  solving  for  each  i  =  0, ...,  m, 


(a  ^ 

v“»J 

where  Uj- 


m 

S  ^kjf  obtaining  the  elements  of  A  as  solutions  of 

k  =  o 


a  ) 

^70 

K^mjJ 

(5.14) 


(5.15) 


for  each  j  =  0,  ...,  n. 

Several  generalizations  and  simplifications  can  be  made  to  this  algorithm.  First 
note  that,  as  stated,  the  matrices  F  and  G  are  square  (i.e.,  have  the  same  number  of  nodes 
as  basis  functions)  to  allow  use  of  the  standard  matrix  inverse  in  solving  Eqs.  5.14  and 
5.15.  The  algorithm  can  be  made  to  work  with  overdetermined  F  and  G  matrices  (more 
nodes  than  basis  functions)  through  use  of  the  pseudoinverse  [Lawson  and  Hanson].  Let 
the  pseudoinverse  of  a  matrix  F  be  denoted  FGL  . 
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Additionally,  the  matrix  A  need  not  be  solved  for  column  by  column.  For  instance, 
one  can  first  define  the  matrix 


^  pi 

S={l(fi,8j)}  ,  (5.16) 

ij  =  0 

where  Mr  and  are  the  number  of  basis  functions  in  each  basis.  Also,  let  the  matrix 

r  Cr 

T 

made  up  of  the  a-  ■  values  be  denoted  as  E.  We  can  then  say  that  GE  =  S  (a  matrix 

(-1)  T 

formulation  of  Eq.  5.14  for  all  /  =  0, ...,  m)  and  then  E  =  G  S  .  Likewise,  we  can 
say  FA  =  E^  (a  matrix  formulation  of  Eq.  5.15  for  all  j  -  0, ...,  n),  and  finally  that 

A  =  (5.17) 

5.3.2.  Use  of  I.P.Q.F.’s  for  the  Extrapolation  Algorithm 

An  I.P.Q.F.  can  be  used  to  approximate  the  integral  involved  in  the  field 
extrapolation  algorithm.  This  is  accomplished  by  using  the  mode  functions  obtained 
through  a  numerically  based  normal  mode  decomposition  algorithm  (such  as  KRAKEN) 
as  the  basis  functions  described  in  the  previous  section  (note  that  since  the  mode  functions 
form  an  orthonormal  set  it  can  also  be  shown  that  they  are  linearly  independent  [Naylor 
and  Sell]).  In  order  for  the  F  and  G  matrices  seen  above  to  be  square  or  overdetermined, 
we  see  that  the  requirement  for  the  number  of  reference  transmitters  is  n  >  M,  where  M  is 
the  number  of  propagating  modes  in  the  problem.  Examples  in  the  following  section  will 
illustrate  the  improvements  in  performance  obtained  using  the  I.P.Q.F. 

By  using  this  integration  method,  we  have  moved  away  from  the  general 
formulation  explored  in  Ch.  2  where  the  modeled  field  p  could  be  generated  through  any 
number  of  acoustic  field  simulation  methods.  While  this  is  still  the  case,  since  the  modal 
decomposition  would  need  to  be  calculated  in  any  case  as  part  of  the  I.P.Q.F.,  it  would 
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make  sense  to  also  use  the  modal  decomposition  to  form  the  modeled  field.  In  this  form 
the  Huygens’  principle-based  formulation  can  be  made  equivalent  to  an  explicit  modal 
decomposition-based  method  of  field  extrapolation  which  was  introduced  in  [Rouseff  et 
al.].  The  modal  decomposition  version  of  the  algorithm  is  an  adaptation  of  processing 
methods  developed  in  [Yang  (1987)]  and  [Yang  (1990)]. 

As  we  did  earlier  in  this  dissertation,  we  will  make  a  simplifying  assumption  that, 
for  cases  where  there  is  a  sedimentary  bottom  layer  into  which  transmit/receive  arrays 
cannot  extend,  the  bottom  parameters  will  be  such  that  any  modes  that  would  propagate 
significant  energy  in  the  bottom  will  be  stripped  off.  In  other  words,  we  will  restrict 
ourselves  to  primarily  water-borne  energy.  This  will  allow  us  to  assume  a  constant  density 
of  p(z)  =  1  for  the  water  column  (a  common  assumption  in  many  cases),  and  to  rewrite 
our  orthonormality  condition  for  the  mode  functions  as 


D 

J  ^Jz)\\f^(z)dz  =  { 
0 


\,m  =  n 
0,m^n 


(5.18) 


where  D  is  the  water  depth. 

5.3.3.  Extrapolation  Examples  Using  I.P.Q.F.’s 

Example  5A  will  use  the  same  basic  scenario  as  Ex.  3E  in  Sec.  3.5,  i.e.,  an 
attenuating  sediment  layer  underlying  the  water  column.  The  only  difference  will  be  that 
instead  of  spacing  transmitters  along  the  reference  array  at  Az  =  1  m,  18  evenly  spaced 
transmitters  will  be  placed  along  the  array  (water  depth  =  150  m,  Az  =  7.89  m )  in  keeping 
with  the  information  displayed  in  Fig.  3.11  which  implied  roughly  18  significantly 
propagating  modes  in  the  problem.  Figure  5.5  shows  the  results  of  using  the  trapezoidal 
integration  method  in  the  extrapolation  algorithm.  The  matching  coefficient  here  is 
approximately  0.956,  which  shows  a  decrease  from  the  1  m  spacing  case  (0.996)  in 
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Figure  5.5.  Ex.  5 A:  trapezoidal  integration. 
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Figure  5.6  shows  the  results  of  using  an  I.P.Q.F.  for  this  scenario.  The  matching 
coefficient  here  is  approximately  1.0  (the  actual  value  calculated  out  to  six  decimal  places 
was  0.999999).  As  can  be  seen  from  the  transmission  loss  comparison  and  the  Argand 
diagram  there  is  virtually  no  difference  between  the  extrapolated  and  ground  truth  fields. 
For  this  I.P.Q.F.  the  matrix  S  of  Eq.  5.16  is  an  18x18  identity  matrix  due  to  the 
orthogonality  of  the  mode  functions.  Since  Eq.  5.12  is  merely  a  better  numerical 
integration  method,  the  obliquity  factors  still  affect  the  result.  Figure  5.6  shows  that  the 
obliquity  factors  (the  obliquity  characterization  number  here  is  approximately  1 .0)  have 
practically  no  effect  for  this  scenario.  The  slight  decrease  in  matching  coefficient  seen  in 
Ex.  3C  can  be  attributed  to  effects  due  to  the  trapezoidal  integration  method  employed 
there. 

Example  5B  uses  the  scenario  of  Ex.  3A,  namely  a  rigid  bottom  with  constant 
sound  speed  in  the  water  column  (1500  m/s).  Since  there  are  10  propagating  modes  for 
this  scenario,  we  will  place  10  equally  spaced  transmitters  along  the  reference  transmit 
array  (water  depth  =  50  m,  Az  =  4.55  m ).  Figure  5.7  shows  the  results  of  using  trapezoidal 
integration  to  do  the  extrapolation.  The  matching  coefficient  here  is  approximately  0.889. 
The  degradation  from  Ex.  3A  is  again  due  to  the  aliasing  effects  due  to  a  lower  sampling 
“rate”  that  are  inherent  in  the  trapezoidal  integration  method. 

Eigure  5.8  shows  the  results  using  an  I.P.Q.F.  By  using  this  method,  the  matching 
coefficient  has  been  raised  to  0.996.  This  is  an  improvement  over  the  value  of  0.982 
obtained  in  Ex.  3A  which  used  trapezoidal  integration  and  transmitters  spaced  1  m  apart. 
We  can  now  say  that  any  degradation  in  the  extrapolated  field  is  due  solely  to  the  obliquity 
factors. 

5.4.  Obliquity  Characterization  Number  vs.  Matching  Coefficient  Revisited 

In  Sec.  3.4  we  derived  an  expression  called  the  obliquity  characterization  number, 
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Table  5.1.  Matching  Coefficients  V5.  Obliquity  Characterization  Numbers. 


Example  # 

Obliquity 
Characterization 
Number  (T ) 

Matching 

Coefficient 

5C  (3C) 

0.999 

0.998 

5B  (3A) 

0.997 

0.996 

5D  (3D) 

0.994 

0.991 

5E  (3E) 

0.983 

0.974 

5F  (3B) 

0.965 

0.952 

Yg,  which  was  intended  as  an  approximation  to  the  matching  coefficient  that  would  be 
seen  in  a  given  environmental  scenario  given  approximations  of  the  obliquity  factors.  We 
saw  there  that  although  there  appeared  to  be  a  linear  relationship  between  this  value  and 
the  matching  coefficient,  there  was  an  offset  such  that  when  the  obliquity  characterization 
number  was  1 .0  (no  obliquity  factor  effects),  the  matching  coefficient  was  somewhat  less 
than  1.0  (-0.988).  From  the  results  of  Sec.  5.3,  we  can  now  see  that  this  was  due  to 
aliasing  in  the  sampling  of  the  field  values  along  the  reference  array.  For  example,  when 
we  use  an  I.P.Q.F.  to  evaluate  the  integral  for  Ex.  3A  (as  in  Ex.  5B),  which  avoids  the 
aliasing  effects,  we  achieve  a  matching  coefficient  of  0.996  vs.  the  value  of  0.982  obtained 
by  using  trapezoidal  integration  with  reference  transmitters  spaced  at  1  m  intervals. 

Table  5.1  shows  the  results  of  using  an  I.P.Q.F.  to  evaluate  the  integrals  for  the 
same  scenarios  that  were  seen  in  Sec.  3.4.  Figure  5.9  shows  a  plot  of  these  new  values 
with  a  least-squares  linear  fit  compared  to  the  trapezoidal  integration  case  (the  sediment 
example  has  been  added  for  this  case  since  the  I.P.Q.F.  has  negated  any  effects  due  to 
element  spacing).  The  equation  for  the  least  squares  fit  of  the  new  data  is 
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Figure  5.9.  Matching  Coejf.  V5.  Obliquity  Char.  Number  and  Linear  Fit. 
element  spacing).  The  equation  for  the  least  squares  fit  of  the  new  data  is 

|^:|^  =1.371X0-0.371.  (5.19) 

We  can  see  that  the  fit  does  extrapolate  to  a  matching  coefficient  of  1.0  for  an  obliquity 
characterization  number  or  1.0.  The  slope  has  become  slightly  steeper,  although  this  can 
most  likely  be  attributed  to  the  clustering  of  data  points  near  the  higher  obliquity 
characterization  numbers.  In  any  case,  the  obliquity  characterization  does  provide  a  useful 
approximation  to  the  expected  matching  coefficient. 

5.5.  Obliquity  Factor  Cancellation  Using  an  I.P.Q.F. 

Our  use  of  the  modal  decomposition  in  evaluating  the  integral  also  allows  us  to 
perfectly  cancel  the  obliquity  factor  effects.  In  Eq.  2.8  we  saw  how  calculating  the  integral 
produced  an  extra  in  the  modal  summation  and  an  incorrect  multiplier  in  front  of 

the  summation.  An  attempt  to  cancel  both  these  effects  was  made  by  multiplying  the 
summation  by  the  mode-independent  factor  N^'  (see  Eq.  2.10).  The  reference 
wavenumber  multiplier  combined  with  the  extra  to  produce  the  modal-dependent 

obliquity  factors.  If  we  modify  our  I.P.Q.F.  by  making 


128 


5  = 


0  0  0  j 

0  ^2^^  0  0 
0  0  ...  0 
0  0  0 


(5.20) 


i.e.,  correct  for  the  extra  factors  that  occur  in  the  integration  by  placing  them  along 

the  diagonal  of  the  5  matrix,  and  subsequently  modify  the  mode-independent  multiplier 
so  that 


''J 


1/2 


(5.21) 


the  effects  of  the  obliquity  factors  should  be  cancelled  exactly.  Note  that  only  the  primary 
extrapolation  now  need  be  done  since  the  alternate  extrapolation  formula  was  developed 
in  order  to  lessen  the  effects  of  the  obliquity  factors.  Figure  5.10  shows  the  results  of 
calculating  the  extrapolation  integral  using  a  “deobliquified”  I.P.Q.F.  As  expected,  the 
matching  coefficient  here  is  1 .0. 

5.6.  Comments  on  Modal  Decomposition  Methods 

Section  5.5  showed  how  modal  decomposition  methods  can  be  used  to  do  field 
extrapolation  (see  also  [Rouseff  cr  a/.]).  There  are  several  interesting  questions  about  these 
methods  that  invite  future  work.  One  area  has  to  do  with  how  elements  are  spaced  on  the 
array.  Although  the  sampling  criterion  for  the  modal  decomposition-based  methods  only 
deals  with  the  actual  number  of  transmitters/receivers  that  must  populate  an  array,  does 
this  mean  that  all  the  elements  could  be  tightly  grouped  in  one  small  vertical  segment  of 
the  array?  Preliminary  numerical  studies  have  shown  that  this  is  not  the  case  due  to 
numerical  instabilities  in  calculating  the  A  matrix.  This  is  a  subject  that  deserves  detailed 
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Study. 

Another  area  of  future  work  focuses  on  modeled  field  mismatch  issues.  Much  of 
the  work  in  Ch.  4  can  translate  directly  to  the  results  of  Ch.  5.  One  extra  consideration, 
however,  is  that  while  the  mode  functions  that  are  used  as  an  integration  basis  are  exact  for 
the  modeled  field,  they  will  not  be  exact  for  the  Green’s  function  data.  As  seen  in  Ch.  4, 
for  environmental  scenarios  in  which  high  order  modes  are  stripped  off,  the  perturbations 
in  the  mode  functions  tend  to  be  small,  so  this  mode  basis  mismatch  should  not  contribute 
much  to  the  errors.  Preliminary  numerical  studies  have  borne  this  out.  This  also  is  a  topic 
that  deserves  further  study. 

5.7.  Summary  of  Chapter  5 

Chapter  5  dealt  with  numerical  integration  issues  pertinent  to  the  acoustic  field 
extrapolation  algorithm  due  to  the  necessity  of  having  discretely  placed  transmitters/ 
receivers  along  the  reference  arrays.  It  was  shown  how  simple  summation-type 
(specifically,  trapezoidal)  numerical  integration  techniques  can  degrade  the  results  for 
insufficiently  dense  sampling  along  the  reference  arrays.  A  method  was  suggested  for 
determining  a  priori  the  required  spacing  for  elements  on  the  reference  arrays.  It  was  then 
shown  how  inner  product  quadrature  formulas  can  be  used,  assuming  a  normal  mode 
expansion  of  the  modeled  field  and  Green’s  function  data,  to  alleviate  this  problem.  The 
sampling  requirement  there  was  that  the  reference  arrays  contain  at  least  as  many  elements 
as  the  number  of  propagating  modes.  It  was  shown  how  using  the  mode  functions  as  the 
integration  basis  essentially  ties  together  the  Huygens’  principle-based  and  modal 
decomposition-based  extrapolation  algorithms. 


Chapter  6:  Conclusion 


6.1.  Summary  of  Dissertation 

This  dissertation  has  concentrated  on  analyzing  algorithms  for  extrapolating 
underwater  acoustic  fields.  These  types  of  algorithms  may  be  necessary  for  cases  where  it 
is  desired  to  predict  acoustic  fields  at  some  distance  from  an  acoustic  source,  but  the 
environmental  conditions  over  the  entire  range  that  are  needed  as  input  to  a  modeling 
algorithm  are  unavailable,  or  measured  data  at  all  locations  are  also  unavailable  [Cox  et 
al].  The  extrapolation  algorithms  combine  modeled  and  measured  data  in  order  to  predict 
the  fields. 

Chapter  1  gave  some  background  material  in  the  form  of  an  introduction  to  normal 
mode  theory  and  a  source  localization  algorithm  (HAP)  that  utilized  medium  calibration 
much  as  the  extrapolation  algorithm.  Chapter  2  outlined  Huygens’  principle  and  the 
Fresnel-Kirchhoff  diffraction  formula  which  then  served  as  models  for  the  underwater 
field  extrapolation  algorithm.  It  was  found  that  by  combining  modeled  fields  along  a 
reference  array  with  calibration  data  taken  along  the  array,  useful  extrapolation  results 
could  be  obtained.  The  combination  of  the  information  was  essentially  an  integral  over 
depth  of  their  product.  By  analyzing  the  extrapolated  and  ground  truth  fields  using  normal 
mode  theory,  it  was  seen  that  the  extrapolated  field  contained  a  modal-dependent 
“obliquity”  factor  which  became  more  prominent  for  more  steeply  propagating  modes.  It 
was  shown  how  an  alternate  extrapolation  formula  based  on  the  range  derivative  of  the 
modeled  field  could  help  average  out  the  obliquity  factors  to  first  order.  It  was  also 
analytically  shown  how  calibration  data  take  into  effect  range  dependencies  in  the 
medium  to  produce  reliable  extrapolated  fields. 

Chapter  3  contained  numerically  modeled  illustrations  of  the  mathematical 
relationships  shown  in  Ch.  2.  Obliquity  factor  effects  were  examined  for  cases  of 
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moderately  steeply,  very  steeply,  and  mostly  horizontally  propagating  energy.  These 
examples  showed  how  more  steeply  propagating  energy  tended  to  degrade  the 
extrapolated  field  estimates  through  a  quantitative  goodness  measure,  the  matching 
coefficient.  It  was  also  shown  how  an  estimate  of  the  degradation  in  matching  coefficient 
could  be  obtained  through  estimates  of  the  obliquity  factors  alone. 

Chapter  4  dealt  with  modeled  field  mismatch  issues,  i.e.,  cases  where  the  assumed 
environmental  parameters  used  as  input  to  the  modeling  step  in  the  algorithm  did  not 
correctly  reflect  the  actual  parameters.  Perturbational  results  were  shown  relating  errors  in 
the  sound  speed  profile  of  the  water  column  to  errors  in  the  horizontal  wavenumbers  and 
mode  functions.  It  was  shown  how  these  errors  translated  into  errors  in  the  extrapolated 
fields.  Perturbational  results  were  also  examined  for  sediment  property  errors  (sound 
speed,  attenuation,  and  density).  As  in  the  obliquity  factor  effects,  it  was  shown  how 
estimates  of  the  eigenvalue  perturbations  alone  can  lead  to  an  estimate  of  the  degradation 
in  matching  coefficient. 

Chapter  5  concentrated  on  the  issue  of  numerically  estimating  the  integral  in  the 
extrapolation  algorithm.  This  was  necessary  due  to  the  fact  that  elements  must  be  placed  at 
discrete  locations  along  the  reference  arrays.  It  was  shown  how  undersampling  the 
modeled  field/calibration  data  product  can  degrade  the  extrapolated  fields  when  simple 
summation-type  (trapezoidal)  quadrature  methods  are  used,  and  how  to  estimate  the 
spacing  required  for  reliable  extrapolations.  It  was  shown  how  inner  product  quadrature 
formulas  could  be  used  to  essentially  perform  a  modal  decomposition-based  extrapolation 
formula  by  using  the  mode  functions  of  the  modeled  field  as  the  integration  basis.  The 
sampling  requirement  was  that  there  be  at  least  as  many  elements  on  the  reference  array  as 
there  were  propagating  modes.  It  was  then  shown  how  the  I.P.Q.F.  could  be  modified  to 
cancel  the  obliquity  factors  exactly. 
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6.2.  Suggestions  for  Future  Work 

There  are  several  areas  for  further  research  in  the  topic  of  field  extrapolation 
algorithms.  One  has  to  do  with  the  effects  of  surface  and  bottom  roughness  on  the 
algorithm.  Bottom  roughness  can  be  difficult  to  model,  but  its  effects  over  the  calibration 
distance  should  be  taken  into  account  by  the  Green’s  function  data,  leaving  only  the 
effects  due  to  bottom  roughness  over  the  extrapolation  range.  Surface  roughness,  unlike 
the  majority  of  bottom  roughness  scenarios,  is  time-variable.  Therefore,  an  area  of  study 
exists  as  far  as  determining  how  calibration  data  are  affected  by  the  changing  sea  surface, 
as  well  as  the  effects  of  coupling  between  surface  and  bottom  roughness.  Along  with  this 
time-variable  theme  is  the  issue  of  internal  waves,  and  how  they  affect  the  calibration  data. 

Another  area  of  study  has  to  do  with  the  out-of-plane  problem.  All  the  material  in 
this  dissertation  assumed  that  the  reference  arrays,  the  hypothetical  source,  and/or  the 
hypothetical  receiver  were  in  a  common  vertical  plane.  The  case  where  they  are  not 
coplanar  deserves  some  attention.  The  problem  wherein  the  reference  arrays  are  not 
perfectly  vertical  provides  another  topic  of  continued  study. 

The  case  of  wideband  signals  is  also  an  area  for  further  research.  This  dissertation 
considered  only  narrowband  cw  signals.  There  may  be  some  advantage  or  disadvantage  to 
using  wideband  signals  for  calibration  of  the  medium. 

Finally,  there  are  the  issues  mentioned  in  Sec.  5.6  having  to  do  specifically  with 
modal  decomposition-based  methods,  namely  the  element  grouping  problem  and  the 
integration  basis  problem  for  mismatch  scenarios. 


Bibliography 


Al-Kurd,  A.  A.,  “Holographic  array  processing  in  the  ocean,”  Ph.D.  dissertation, 
University  of  Washington,  Seattle,  WA,  1993.  Also:  same  author  and  title,  APL-UW  TR 
9306,  Applied  Physics  Laboratory,  University  of  Washington,  Seattle,  1993. 

Baggeroer,  A.  B.,  W.  A.  Kuperman,  and  P.  N.  Mikhalevsky,  “An  overview  of  matched  field 
methods  in  ocean  acoustics,”  IEEE  J.  Oceanic  Eng.,  vol.  18,  no.  4,  pp.  401-424,  1993. 

Boyles,  C.  Py.,  Acoustic  Waveguides  -  Applications  to  Oceanic  Science,  New  York:  John 
Wiley  &  Sons,  1984. 

Brekhovskikh,  L.  M.  and  Y.  P.  Lysanov,  Fundamentals  of  Ocean  Acoustics.  New  York: 
Springer- Verlag,  1991. 

Collins,  M.  D.,  “FEPE  User’s  Guide,”  NORDA  Technical  Note  365,  Naval  Research 
Laboratory,  Stennis  Space  Center,  Mississippi,  October,  1988.  (For  the  FEPE  source  code, 
the  author  may  be  contacted  by  e-mail:  COLLINS@v5160.nrl.navy.mil) 

Collins,  M.  D.  and  E.  K.  Westwood,  “A  higher  order  energy  conserving  parabolic  equation 
for  range  dependent  ocean  depth,  sound  speed,  and  density,”  J.  Acoust.  Soc.  Am.,  vol.  89, 
no.  3,  pp.  1068-1075,  1991. 

Cox,  H.,  H.  Lai  and  M.  Hirano,  “Reciprocity  based  channel  compensation  for  wideband 
communications  in  a  multipath  environment,”  in  Conference  Record  of  the  Twenty- 
Seventh  Asilomar  Conference  on  Signals,  Systems,  and  Computers  (Pacific  Grove,  CA), 
Nov.  1993,  pp.  593-597. 

Davis,  P.  J.  and  P.  Rabinowitz,  Methods  of  Numerical  Integration.  Orlando:  Academic 
Press,  Inc.,  1984. 

Dosso,  S.  E.  and  N.  R.  Chapman,  “Measurement  and  modeling  of  downslope  acoustic 
propagation  loss  over  a  continental  slope,”  J.  Acoust.  Soc.  Am.,  vol.  81,  no.  2,  pp.  258-268, 
1987. 

Dosso,  S.E.,  M.  L.  Yeremy,  J.  M.  Ozard  and  N.  R.  Chapman,  “Estimation  of  ocean-bottom 
properties  by  matched-field  inversion  of  acoustic  field  data,”  IEEE  J.  Oceanic  Eng.,  vol. 
18,  no.  3,pp.  232-239,  1993. 

Dowling,  D.R.,  “Acoustic  pulse  compression  using  passive  phase-conjugate  processing,” 
J.  Acoust.  Soc.  Am.,  vol.  95,  no.  3,  pp.  1450-8,  1994. 


135 


Evans,  R.  B.,  “Coupled  mode  solution  for  acoustic  propagation  in  a  waveguide  with 
stepwise  depth  variations  of  a  penetrable  bottom,”  J.  Acoust.  Soc.  Am.,  vol.  74,  no.  1,  pp. 
188-195,  1983. 

Fox,  W.  L.  J.,  D.  Rouseff  and  R.  R  Porter,  “Underwater  acoustic  field  extrapolation  using 
vertical  reference  arrays,”  in  Conference  Record  of  the  Twenty-Seventh  Asilomar 
Conference  on  Signals,  Systems,  and  Computers  (Pacific  Grove,  CA),  Nov.  1993,  pp.  154- 
158. 

Gingras,  D.F.,  “Impact  of  uncertain  environmental  knowledge  on  the  shallow-water 
transfer  function,”  SACLANTCEN  SM-263,  SACLANT  Undersea  Research  Centre,  La 
Spezia,  Italy,  1992. 

Goncharov,  V.  V.  and  A.  G.  Voronovich,  “An  experiment  on  matched-field  acoustic 
tomography  with  continuous  wave  signals  in  the  Norway  Sea,”  J.  Acoust.  Soc.  Am., 
vol.  93,  no.  4,  pp.  1873-1881,  1993. 

Goodman,  J.  W.,  Introduction  to  Fourier  Optics.  San  Francisco:  McGraw  Hill,  1968. 

Gribble,  J.  D.,  “Interpol atory  inner  product  quadrature  formulas,”  BIT,  vol.  20,  pp.  466- 
474,  1980. 

Halliday,  D.  and  R.  Resnick,  Physics.  New  York:  Wiley,  1978. 

Hamilton,  E.  L.,  “Geoacoustic  modeling  of  the  sea  floor,”  J.  Acoust.  Soc.  Am.,  vol.  68,  no. 
5,pp.  1313-1340,  1980. 

Hamson,  R.  M.  and  R.  M.  Heitmeyer,  “Environmental  and  system  effects  on  source 
localization  in  shallow  water  by  the  matched-field  processing  of  a  vertical  array,” 
J.  Acoust.  Soc.  Am.,  vol.  86,  no.  5,  pp.  1950-1959,  1989. 

Jensen,  F.  B.,  W.  A.  Kuperman,  M.  B.  Porter,  and  H.  Schmidt,  Computational  Ocean 
Acoustics.  New  York:  American  Institute  of  Physics,  1994. 

Kinsler,  L.  E.,  A.  R.  Frey,  A.  B.  Coppens,  and  J.  V.  Sanders,  Fundamentals  of  Acoustics. 
New  York:  John  Wiley  and  Sons,  1982,  3rd  ed. 

Lawson,  C.  L.  and  R.  J.  Hanson,  Solving  Least  Squares  Problems.  Englewood  Cliffs,  NJ: 
Prentice-Hall,  1974. 

Lindsay,  C.  E.  and  N.  R.  Chapman,  “Matched  field  inversion  for  geoacoustic  model 
parameters  using  adaptive  simulated  annealing,”  IEEE  J.  Oceanic  Eng.,  vol.  18,  no.  3,  pp. 
224-231,  1993. 


136 


Marks,  R.  J.  II,  Introduction  to  Shannon  Sampling  and  Interpolation  Theory.  New  York: 
Springer- Verlag,  1991. 

McGrath,  J.  R,  “Gaussian  product-type  quadratures,”  App/.  Math,  and  Comp.,  vol.  5,  pp. 
265-280,  1979. 

Mourad,  R  D.,  D.  Rouseff,  R.  P.  Porter,  and  A.  Al-Kurd,  “Source  localization  using  a 
reference  wave  to  correct  for  oceanic  variability,”  J.  Acoust.  Soc.  Am.,  vol.  92,  no.  2,  pp. 
1031-1039,  1992. 

Naylor,  A.  W.  and  G.  R.  Sell,  Linear  Operator  Theory  in  Engineering  and  Science.  New 
York:  Springer- Verlag,  1982. 

Oppenheim,  A.  V.  and  R.  W.  Schafer,  Discrete-Time  Signal  Processing.  Englewood  Cliffs, 
New  Jersey:  Prentice  Hall,  1989. 

Oppenheim,  A.  V.,  A.  S.  Willsky  and  I.  T.  Young,  Signals  and  Systems.  Englewood  Cliffs, 
New  Jersey:  Prentice  Hall,  1983. 

Papoulis,  A.,  Probability,  Random  Variables,  and  Stochastic  Processes.  New  York: 
McGraw-Hill,  1984. 

Porter,  M.  B.,  “The  KRAKEN  Normal  Mode  Program,”  Unpublished  technical  report, 
S  ACL  ANT  Undersea  Research  Centre,  11  Nov.  1992.  (The  source  code  for  KRAKEN  is 
available  via  anonymous  FTP  from  alba.njit.edu) 

Porter,  R.  P,  “Three-Dimensional  Imaging  Methods  in  Acoustics  and  Electromagnetics,” 
Unpublished  lecture  notes  for  EE  578,  Dept,  of  Electrical  Engineering,  University  of 
Washington,  Seattle,  1994. 

Porter,  M.  B.  and  E.  L.  Reiss,  “A  numerical  method  for  ocean  acoustic  normal  modes,”  J. 
Acoust.  Soc.  Am.,  vol.  76,  no.  1,  pp.  244-252,  1984. 

Porter,  R.  R,  P.  D.  Mourad  and  A.  Al-Kurd,  “Wave-front  reconstruction  in  variable, 
multimode  waveguides,”  J.  Opt.  Soc.  Am.  A,  vol.  9,  no.  9,  pp.  1984-1990,  1992. 

Porter,  R.  P,  D.  Rouseff,  W.  L.  J.  Fox  and  M.  Siderius,  “Acoustic  calibration  in  shallow 
water:  Theory,  simulations,  and  a  preliminary  site  study,”  APL-UW  TM  31-93,  Applied 
Physics  Laboratory,  University  of  Washington,  Seattle,  1993. 

Press,  W.  H.,  B.  P.  Flannery,  S.  A.  Teukolsky,  and  W.  T.  Vetterling,  Numerical  Recipes 
in  C.  New  York:  Cambridge  University  Press,  1988. 


137 


Rajan,  S.  D.,  J.  F.  Lynch  and  G.  V.  Frisk,  “Perturbative  inversion  methods  for  obtaining 
bottom  geoacoustic  parameters  in  shallow  water,”  J.  Acoust.  Soc.  Am.,  vol.  82,  no.  3,  pp. 
998-1017,  1987. 

Rouseff,  D.,  “Ocean  acoustic  holography:  Using  a  reference  source  to  remove  oceanic 
variability,”  Unpublished  technical  report.  Applied  Physics  Laboratory,  University  of 
Washington,  Seattle,  1989. 

Rouseff,  D.,  (Untitled),  Unpublished  technical  note.  Applied  Physics  Laboratory, 
University  of  Washington,  Seattle,  February  8,  1993. 

Rouseff,  D.,  (Untitled),  Unpublished  technical  note.  Applied  Physics  Laboratory, 
University  of  Washington,  Seattle,  April  29,  1993. 

Rouseff,  D.,  M.  Siderius,  W.  L.  J.  Fox  and  R.  P.  Porter,  “Acoustic  calibration  in  shallow 
water  using  sparse  data,”  submitted  and  accepted  for  publication  in  J.  Acoust.  Soc.  Am., 
August,  1994. 

Schiff,  L.  I.,  Quantum  Mechanics.  New  York:  McGraw-Hill,  1968. 

Uscinski,  B.  J.,  The  Elements  of  Wave  Propagation  in  Random  Media.  New  York: 
McGraw-Hill,  1977. 

Wen,  T.,  W.  J.  Felton,  J.  C.  Luby,  W.  L.  J.  Fox  and  K.  L.  Kientz,  “Environmental 
Measurements  in  the  Beaufort  Sea,  Spring  1988,”  APL-UW  TR  8822,  Applied  Physics 
Laboratory,  University  of  Washington,  Seattle,  WA,  March  1989. 

Yang,  T.  C.,  “A  method  of  range  and  depth  estimation  by  modal  decomposition,”  J. 
Acoust.  Soc.  Am.,  vol.  82,  no.  5,  pp.  1736-1745,  1987. 

Yang,  T.  C.,  “Effectiveness  of  mode  filtering:  a  comparison  of  matched-field  and  matched¬ 
mode  processing,”  J.  Acoust.  Soc.  Am.,  vol.  87,  no.  5,  pp.  2072-2084,  1990. 

Yariv,  A.,  An  Introduction  to  Theory  and  Applications  of  Quantum  Mechanics.  New  York: 
John  Wiley  and  Sons,  Inc.,  1982. 


Vita 


Warren  Leonard  John  Fox  was  bom  March  30,  1966,  in  Seattle,  Washington.  He 
grew  up  in  Everett,  Washington,  with  his  father  John,  mother  Karen,  brother  Colin,  and 
sisters  Sarah  and  Anne.  He  graduated  from  Cascade  High  School  in  Everett,  Washington, 
in  June  of  1984.  He  was  co-valedictorian  of  his  graduating  class,  a  National  Merit  Scholar, 
and  the  1984  Northwest  District  high  jump  champion.  In  June  of  1988  he  graduated  cum 
laude  with  his  Bachelor  of  Science  in  Electrical  Engineering  degree  from  the  University  of 
Washington  in  Seattle,  Washington.  Upon  graduation,  he  joined  the  professional  staff  of 
the  Applied  Physics  Laboratory,  University  of  Washington,  as  an  electrical  engineer.  In 
September  of  1989  he  was  granted  an  APL-UW  Fellowship  to  support  his  graduate  work. 
He  earned  his  Master  of  Science  in  Electrical  Engineering  from  the  UW  in  December  of 
1990,  where  his  thesis  topic  was  the  use  of  a  new  time-frequency  representation  for  range- 
Doppler  applications.  He  continued  on  the  APL-UW  Fellowship  during  his  Ph.D.  studies. 
His  research  interests  include:  signal  processing,  detection/classification/localization 
applications,  and  ocean  acoustic  applications.  He  is  a  member  of  the  Signal  Processing 
and  Oceanic  Engineering  Societies  of  the  Institute  of  Electrical  and  Electronic  Engineers 
(I.E.E.E.). 

Selected  Publications: 


Fox,  W.  L.  J.,  “Application  of  a  Cone-Shaped  Kernel  Time-Frequency  Representation  to 
Sonar  and  Radar  Range-Doppler  Processing,”  Master’s  thesis.  University  of  Washington, 
Seattle,  WA,  Dec.  1990. 

Fox,  W.  L.  J.,  D.  Rouseff  and  R.  P.  Porter,  “Underwater  acoustic  field  extrapolation  using 
vertical  reference  arrays,”  in  Conference  Record  of  the  Twenty-Seventh  Asilomar 
Conference  on  Signals,  Systems,  and  Computers  (Pacific  Grove,  CA),  Nov.  1993,  pp.  154- 
158. 


REPORT  DOCUMENTATION  PAGE 


Form  Approv&d 
OPMNo,  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  reviewing  the  collection  of  Information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggestions 
for  reducing  this  burden,  to  Washington  Headquarters  Services.  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway.  Suite  1204,  Arlington,  VA  22202-4302,  and  to 
the  Office  of  Information  and  Regutatory  Affairs,  Office  of  Management  and  Budget.  Washington.  DC  20503. 


2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

March  1995  Technical 


1.  AGENCY  USE  ONLY  (Leave  blank) 


4.  TITLE  AND  SUBTITLE 

Underwater  Acoustic  Field  Extrapolation:  Theory  and  Sensitivity  Analysis 


6.  AUTHOR{S) 

Warren  L.  J.  Fox 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS{ES) 
Applied  Physics  Laboratory 
University  of  Washington 
1013  NE  40th  Street 
SeatUe,  WA  98105-6698 


9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

APL-UW  Independent  Research  & 

Development  Program 
1013  NE40th  Street 
SeatUe,  WA  98105-6698 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

APL-UW  TR  9501 


10.  SPONSORING  /  MONITORING 
AGENCY  REPORT  NUMBER 


11.  SUPPLEMENTARY  NOTES 

12a.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Distribution  Unlimited 

12b.  DISTRIBUTION  CODE 

13.  ABSTRACT  (Maximum  200  words) 

In  some  underwater  acoustic  applications,  it  is  desirable  to  predict  what  an  acoustic  field  will  be  at  some  distance  away 
from  a  source.  Numerical  modeling  base^  on  environmental  parameters  (e.g.,  sound  speed,  bathymetry,  sediment  properties)  is  one 
possibility,  but  the  amount  of  environmental  information  necessary  for  accurate  long-range  modeling  is  often  beyond  what  is  real¬ 
istically  measurable.  Actual  acoustic  measurements  are  another  possibility,  but  our  knowledge  of  the  fields  is  then  limited  to  the 
points  in  the  ocean  where  the  sources  and  receivers  are  located.  We  present  a  hybrid  method  for  predicting  acoustic  fields  that  lakes 
limited  acoustic  measurements  and  extrapolates  them  over  short  ranges  using  modeled  fields.  The  measured  and  modeled  fields  are 
combined  in  an  integral  formulation  via  a  specialization  of  Huygens’  principle.  The  basic  theoretical  formulation  is  stated  and  ana¬ 
lyzed  using  the  theory  of  normal  modes.  The  formulation  leads  directly  to  modal-dependent  (grazing  angle  dependent)  “obliquity” 
factors,  which  are  canceled  to  first  order  by  averaging  with  a  range  derivative-based  extrapolation.  Sensitivity  of  the  algorithm  to 
the  higher  order  components  of  the  obliquity  factors  is  studied.  As  with  other  model-based  algorithms,  environmental  mismatch  can 
degrade  the  results.  The  algorithm’s  sensitivity  to  mismatch  in  water  column  sound  speed  and  sediment  sound  speed,  attenuation, 
and  density  are  studied  via  a  perturbaiional  approach.  Since  the  algorithm  explicitly  contains  an  integral  and  the  measured  field  is 
collected  at  discrete  locations,  numerical  quadrature  techniques  are  studied.  It  is  shown  how  undersampling  can  degrade  the  qual¬ 
ity  of  the  extrapolated  fields,  and  how  the  spatial  sampling  required  for  reliable  extrapolation  can  be  obtained.  A  modal  decompo¬ 
sition-based  integration  method  is  then  shown,  where  the  number  of  reference  elements  required  is  equal  to  the  number  of  modes 
that  contribute  significantly  to  the  fields.  A  modified  version  of  this  formulation  can  exactly  cancel  the  obliquity  factors. 


14.  SUBJECT  TERMS 

Underwater  Acoustics,  Acoustic  Modeling,  Normal  Mode  Theory,  Huygens’  Principle, 
Obliquity  Factors,  Environmental  Mismatch,  Perturbation  Theory,  Inner  Product  Quadrature 
Formulas 


17.  SECURITY  CLASSIFICATION 
OF  REPORT 

Unclassified 


18.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 

Unclassified 


19.  SECURITY  CLASSIFICATION 
OF  ABSTRACT 

Unclassified 


15.  NUMBER  OF  PAGES 

149 


16.  PRICE  CODE 


20.  LIMITATION  OF  ABSTRACT 


NSN  7540-01-280-5500 


Standard  Form  298  (Rev.  2-89) 
Prescribed  by  ANSI  Std.  239-18 
299-01 


